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I.  INTRODUCTION 

This  report  presents  a  satellite  simulation  algorithm  written  for  the  IBM  PC  and  clones 
in  Microsoft  QuickBASIC  4.5.  The  program  can  either  generate  an  ephemeris  of  satellite 
ground-track  position  and  location  with  respect  to  a  fixed  ground  station,  or  determine  the 
position  and  time  the  satellite  is  above  the  horizon  of  a  specified  ground  station.  Provision 
is  made  for  the  program  to  handle  multiple  satellites  and  multiple  ground  stations.  Input 
is  either  by  way  of  the  computer  keyboard  or  by  separate  files  for  satellites  and  ground 
stations.  Output  is  written  both  on  the  screen  and  to  an  ASCII  output  file  which  can  be 
accessed  by  a  data  based  management  program. 

Two  versions  of  the  satellite  kinematics  subroutine  pos. update  are  presented.  These 
algorithms  were  provided  by  NAVSPASUR  [Ref.  1].  The  first  is  a  very  accurate  algorithm 
based  upon  a  paper  by  Brouwer  [Ref.  2]  and  is  contained  in  NAVSPASUR's  FORTRAN  pro- 
gram called  SHOWALL.  A  listing  of  the  BASIC  version  of  this  algorithm  is  presented  in 
Appendix  E.  The  second  algorithm  is  a  much  simplified  approximation  to  the  Brouwer 
model.  Equations  for  the  simplified  model  are  contained  in  Appendix  A,  and  the  BASIC 
implemention  is  embedded  in  the  listing  of  the  full  simulation  program  in  Appendix  D. 

The  next  sections  contain  program  operating  instructions,  structure  of  the  input/ 
output  files,  and  suggestions  for  future  improvements.  The  computation  of  the  times  of 
culmination,  rise  and  set  are  in  Appendix  B,  and  the  computation  of  satellite  heading  is  in 
Appendix  C.  Readers  unfamiliar  with  the  nomenclature  of  satellite  orbital  elements  should 
consult  any  standard  work  such  as  Escobal  [Ref.  3,  Ch.  3]. 

II.  OPERATLNG  INSTRUCTIONS 

A.  Overview.  The  source  code  is  named  satsta.bas  and  the  executable  code  is  named 
satsta.exe.  To  run,  enter  SATSTA  at  the  DOS  prompt  and  press  the  <_i  key. 

The  first  menu  (Fig.  1)  requests  a  data  input  mode.  Pressing  1  or  2  will  require  the 
data  for  one  satellite  and  one  ground  station  to  be  input  from  the  keyboard.  Pressing  3  or  4 
assumes  that  two  data  files  are  accessible:  either  elements.dat  or  elements. nss  containing 
the  orbital  elements  for  one  or  more  satellites,  and  station.dat  containing  the  location  of 
one  or  more  ground  stations — the  structure  of  these  files  is  discussed  in  the  next  section. 
No  matter  which  of  these  four  options  is  chosen,  a  simulation  starting  time,  ending  time 
and  time-step  increment  will  be  required  input  from  the  keyboard  in  the  third  menu. 

The  second  menu  (Fig.  2)  requests  a  run  mode.  Option  1  creates  an  ephemeris  of 
satellite  geographic  position  and  satellite  position  with  respect  to  the  ground  station  for 
equally  spaced  times  for  the  entire  simulation  duration.    Option  2  computes  the  satellite 
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geographic  position  and  satellite  position  with  respect  to  the  ground  station  only  for  those 
times  during  the  simulation  duration  at  which  the  satellite  is  above  the  horizon  of  the 
ground  station. 

The  third  menu  (Fig.  3)  requests  the  simulation-times,  consisting  of  a  starting  time, 
an  ending  time,  and  a  time-step  increment  (At). 

The  fourth  and  fifth  menus  (Figs.  4  and  5,  respectively)  appear  only  if  the  first  menu 
(data  input  mode)  options  1  or  2  were  chosen.  Menu  four  requests  the  satellite  epoch  and 
elements,  and  menu  five  requests  the  station  coordinates. 

No  matter  which  options  are  chosen,  output  is  directed  to  the  screen  and  is  appended 
to  file  output.dat.  If  the  file  does  not  exist,  it  is  created.  This  ASCII  file  is  structured  so  that 
it  can  be  used  as  input  to  a  data  base  management  program,  if  desired.  The  structure  of 
this  file  is  discussed  in  a  later  section. 

B.  Input.  Menu  1  requests  the  selection  of  a  data  input  mode.  Depending  upon  the 
source  of  satellite  elements,  there  are  two  ways  in  which  these  data  can  be  input.  The 
first  is  using  standard  "epoch  of  date"  elements  in  which  the  time  of  perigee  passage  is  the 
epoch.  The  second  is  the  NAVSPASUR  "One-Line  Charlie"  elements  in  which  the  longitude 
(or  right  ascension)  of  the  ascending  node  has  been  referenced  to  0000  hours,  1  January 
19S5.  The  "One-Line  Charlie"  option  has  been  added  for  those  users  who  may  obtain 
satellite  elements  from  NAVSPASUR.  An  example  of  "epoch  of  date"  input  is  shown  in 
Fig.  3. 

Select  data  input  mode  1,  2,  3  or  4: 

Keyboard  input  for  satellite  ft  ground  station: 

1 .  Epoch  of  date . 

2.  NAVSPASUR  Reference  Date. 

Data  file  input  for  satellite  k   ground  station: 

3.  Epoch  of  date. 

4.  NAVSPASUR  One-Line  Charlie. 

Figure  1.  Menu  1 — Select  data  input  mode. 

If  data  file  input  Option  3  is  chosen  from  Menu  1,  two  ASCII  data  files,  elements  .DAT  and 
station.dat,  must  have  been  previously  created.  The  elements.dat  file  must  contain  a  single 
line  of  for  each  satellite  of  interest.  Each  line  must  contain  10  entries,  each  enclosed  in 
double  quotations  and  separated  by  a  blank  space  (Fig.  6).  The  entries  are:  (1)  a  satellite 
identification  consisting  of  at  most  five  alphanumeric  characters;  (2)  the  epoch  date  (date 
of  perigee  passage)  in  the  form  yyyy  .mmdd,  where  yyyy  is  the  year,  mm  is  the  month,  and 
dd  is  the  day;  (3)  the  epoch  time  (time  of  perigee  passage)  in  the  form  hh.mmss,  where  hh 
is  the  hour,  mm  is  the  minutes,  and  ss  is  the  seconds;  (4)  the  orbital  eccentricity;  (5)  the 


longitude  (or  right  ascension)  of  the  ascending  node  in  degrees;  (6)  the  orbital  inclination 
in  degrees:    (7)  the  argument  of  perigee  in  degrees;    (8)  the  mean  anomaly  in  degi'  • 
(9)  the  mean  motion  in  revolutions  per  day;  and  (10)  the  orbital  decay  in  revolutions  per 
day- squared. 

If  data  file  input  Option  4  is  chosen  from  Menu  1,  two  ASCII  data  files,  elements. nss 
and  station.dat,  must  have  been  previously  created.  The  elements. nss  file  must  contain  a 
single  line  of  for  each  satellite  of  interest  in  the  NAVSPASUR  format.  Each  line  must  contain 
10  entries  starting  in  column  1  and  ending  in  column  65  (Fig.  7).  The  entries  are:  (1)  a 
satellite  identification  consisting  of  at  most  five  alphanumeric  characters  (columns  1-5); 
(2)  the  mean  anomaly  in  revolutions  (columns  6-13);  (3)  the  mean  motion  in  radians  per 
Herg1  (columns  14-21);  (4)  the  orbital  decay  in  radians  per  Herg-squared  (columns  22- 
27);  (5)  the  orbital  eccentricity  (columns  28-35);  (6)  the  argument  of  perigee  in  revolutions 
(columns  36-43);  (7)  the  longitude  (or  right  ascension)  of  the  ascending  node  in  revolutions 
(columns  44-51):  (S)  the  orbital  inclination  in  revolutions  (columns  52-59);  and  (9)  the 
epoch  date  (date  of  perigee  passage)  in  the  form  yymmdd,  where  yyyy  is  the  year,  mm  is  the 
month,  and  dd  is  the  day  (columns  60-65).  A  decimal  point  must  be  present  in  columns  G. 
14,  22,  28,  3G.  44  and  52. 

Select    run    mode    1    or    2: 

1.  Create   an   ephemeris. 

2.  Find   times    satellite   is    above   the   horizon. 

Figure  2.   Menu  2 — Select  run  mode. 

The  station.dat  file  must  contain  a  single  line  of  for  each  ground  station  or  geographical 
location  of  interest.  Each  line  must  contain  5  entries,  each  enclosed  in  double  quotations 
and  separated  by  a  blank  space  (Fig.  S).  The  entries  are:  (1)  a  station  identification 
consisting  of  at  most  three  alphanumeric  characters;  (2)  the  station  latitude  in  the  form 
dd.mmss  (minus  if  south),  where  dd  is  degrees,  mm  is  minutes,  and  ss  is  seconds;  (3)  the  sta- 
tion longitude  in  the  form  ddd.mmss  (minus  if  west);  (4)  the  station  altitude  in  feet  (minus 
if  below  sea  level);  and  (5)  the  name  of  the  station  or  station  location  (alphanumeric). 

If  a  data  file  input  option  is  chosen,  then  SATSTA  will  open  the  two  files  described 
above  and  compute  either  an  ephemeris  or  the  times  that  each  satellite  is  above  the  horizon 
(Option  1  or  2  of  the  second  menu)  for  each  ground  station.  Otherwise,  keyboard  input  for 
one  satellite  and  one  ground  station  will  be  requested  from  the  fourth  and  fifth  menues, 
respectively. 


The  Herg— a  unit  of  time  named  after  the  astronomer  Paul  Herget — is  defined  as  the  period,  divided  by 
27T,  of  a  fictitious  satellite  at  zero  altitude  over  a  smooth  spherical  earth  with  no  atmosphere. 


Simulation-Time   Data: 

Starting   date    (yyyy .nundd)  ?    1983.0201 

Starting   time    (hh.mmss)  ?   00.00 

Ending  date    (yyyy.mmdd)  ?   1983.0202 

Ending   time    (hh.mmss)  ?   00.00 

Time   increment    (minutes)  ?  01 

Press   any  key  to   continue. 
Figure  3.  Menu  3 — Select  simulation  times. 

Five  pieces  of  information  regarding  the  simulation  are  requested  from  the  third  menu 
(Fig.  3):  (1)  the  simulation  starting  date  in  the  form  yyyy.mmdd,  (2)  the  simulation  starting 
time  in  the  form  hh.mmss;  (3)  the  simulation  ending  date  in  the  form  yyyy.mmdd;  (4)  the 
simulation  ending  time  in  the  form  hh.mmss;  and  (5)  the  simulation  time-step  increment 
in  minutes.  The  screen  will  display  sample  data  for  each  entry  in  the  appropriate  format. 
If  any  sample  entry  is  to  be  changed,  it  should  be  completely  overwritten  in  order  for  it 
to  be  read  in  properly.  When  each  entry  is  correct,  enter  it  into  the  computer  by  pressing 
the  «_J  key. 

If  the  keyboard  entry  option  was  chosen  from  the  first  menu,  then  two  more  menues 
will  appear,  otherwise  the  simulation  will  commence  with  input  being  obtained  from  the 
two  data  files. 

Satellite  Data: 

Epoch  date  (yyyy.mmdd)      ?  1983.0201 
Epoch  time  (hh.mmss)        ?  00.00 

Input  mean  elements  of  epoch: 

Satellite  ID  Number  ?  11111 

Eccentricity  ?  0.0005545 

Long,  ascending  node  (deg)   ?  272.43497 

Inclination  (deg)  ?  65.06057 

Arg.  of  perigee  (deg)  ?  295.41470 

Mean  anomaly  (deg)  ?  258.10682 

Mean  motion  (rev/day)  ?  15.44194 

Decay  (rev/day"2)  ?  0 

Press  any  key  to  continue. 

Figure  4.  Menu  4 — Input  satellite  data. 

Menu  4  (Fig.  4)  requires  ten  separate  entries:  (1)  the  epoch  date  (date  of  perigee 
passage)  in  the  form  yyyy.mmdd,  where  yyyy  is  the  year,  mm  is  the  month,  and  dd  is  the 
day;  (2)  the  epoch  time  (time  of  perigee  passage)  in  the  form  hh.mmss,  where  hh  is  the 
hour,  mm  is  the  minutes,  and  ss  is  the  seconds;  (3)  a  satellite  identification  consisting  of  at 
most  five  alphanumeric  characters;  (4)  the  orbital  eccentricity;  (5)  the  longitude  (or  right 
ascension)  of  the  ascending  node  in  degrees;  (6)  the  orbital  inclination  in  degrees;  (7)  the 
argument  of  perigee  in  degrees;  (8)  the  mean  anomaly  in  degrees;  (9)  the  mean  motion  in 


revolutions  per  day;  and  (10)  the  orbital  decay  in  revolutions  per  day-squared.  The  screen 
will  display  sample  data  for  each  entry  in  the  appropriate  format.  Again,  if  any  sample 
entry  is  to  be  changed,  it  should  be  completely  overwritten  in  order  for  it  to  be  read  in 
properly.  When  each  entry  is  correct,  enter  it  into  the  computer  by  pressing  the  <_i  key. 

Station   Coordinates: 

Station   Identification  Number  7  001 

Station  Location   or  Name  ?  Daisy,    Tenn. 

Station   latitude      (dd.mmss,    South  minus)    ?  35.12 

Station   longitude    (ddd.mmss,    West  minus)    ?  -85.12 

Station   altitude      (feet)  ?  500.0 

Press   any   key   to   continue. 

Figure  5.   Menu  5 — Input  ground  station  data. 

Menu  5  (Fig.  5)  requires  five  separate  entries:  (1)  a  station  identification  consisting  •: 
at  most  three  alphanumeric  characters;  (2)  a  station  location  or  name:  (3)  the  station  lati- 
tude in  the  form  dd.mmss  (minus  if  south);  (4)  the  station  longitude  in  the  form  ddd.mmss 
(minus  if  west):  and  (5)  the  station  altitude  in  feet  (minus  if  below  sea  level). 

C.  Output.  Output  is  written  to  the  screen  and  is  appended  to  the  data  file  output.dat 
if  it  exists,  otherwise  file  output.dat  is  created.  A  sample  from  the  output.dat  file  is  shown 
in  Fig.  0 — the  structure  of  this  file  is  discussed  later. 

Sample  output  from  the  "create  an  ephemeris"  option  screen  is  shown  in  Fig.  10. 
The  apogee/perigee  output  is  approximately  valid  only  for  the  first  orbit.  In  particular. 
the  longitude  is  not  corrected  for  the  earth's  rotation  during  the  current  orbit.  Also,  the 
ascending  node  and  argument  of  perigee  are  not  corrected  for  their  rates  of  change  during 
the  current  orbit.  It  should  be  further  noted  that  for  nearly  circular  orbits,  the  height  of 
the  satellite  will  probably  get  lower  than  the  perigee  value  and  higher  than  the  apogee 
value  at  locations  other  than  the  perigee  and  apogee  positions — this  is  an  effect  of  the 
earth's  oblateness. 

The  main  body  of  the  ephemeris  output  consists  of  date,  time,  satellite  latitude, 
longitude  and  height,  followed  by  the  elevation  angle,  azimuth,  and  range  (in  kilometers) 
of  the  satellite  from  the  station.  The  satellite  "look  angle'1  is  the  angle  between  the  satellite 
"ground-zero"  point  (usually  called  the  satellite  "sub-point")  and  the  station.  The  last 
column  of  output  is  the  satellite  heading. 

A  sample  of  the  output  from  the  "times  satellite  is  above  the  horizon"  option  screen 
for  a  time-step  (At)  of  1  minute  is  shown  in  Fig.  11.    First  the  time  of  culmination  for 
the  current  orbit  is  computed  and  saved.   The  culmination  time  is  then  incremented  i 
decremented  bv  At  until  the  satellite  is  below  the  horizon.    The  time  of  rise  and  set  is 
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Satellite:  11111 


Station:  001    Daisy,  Tenn 


Apogee : 
Height 
Latitude 
Longitude 

Perigee : 
Height 
Latitude 
Longitude 


442.8  km. 
-55.1540 

100.2099 


450.4   km. 

55.1539 
280.2099 


Tue  1  Feb  1983 


height 

look 

year 

mo  da  hr 

mn 

sec 

lat 

long 

(km) 

el 

JV 

azi.ii 

range 

angle 

head 

1983 

2   . 

L   0 

0 

0.0 

-12.26 

327.56 

434.0 

-31 

60 

123.30 

7437 

52.90 

157.82 

1983 

2  : 

L   0 

10 

0.0 

-45.84 

347.68 

441.5 

-50 

59 

136.39 

10382 

36.56 

145.02 

1983 

2  : 

L   0 

20 

0.0 

-65.19 

48.39 

446.4 

-69 

07 

151.85 

12354 

19.68 

88.46 

1983 

2  : 

L   0 

30 

0.0 

-44.71 

106.99 

439.1 

-83 

25 

221.39 

13089 

6.45 

34.07 

1983 

2  : 

.   0 

40 

0.0 

-10.95 

126.55 

430.7 

-70 

51 

302.57 

12481 

18.11 

22.04 

1983 

2  : 

L     0 

50 

0.0 

24.08 

141.07 

435.9 

-51 

76 

316.67 

10574 

35.26 

24.38 

1983 

2  : 

L   1 

0 

0.0 

55.55 

168.97 

449.2 

-32 

06 

324.49 

7561 

52.31 

46.47 

Figure  10.  Output  from  ephemens  option.  Epoch  of  date. 


Satellite:  11111 


Station:  001    Daisy,  Tenn. 


year  mo  da  hr  mn   sec    lat    long 


height 
(km)    elev 


look 
azim  range  angle   head 


1983 
1983 
1983 
1983 
1983 
1983 
1983 
1983 
1983 
1983 
1983 
1983 
1983 


1  13 
1  13 


14 

15 
16 
17 


19.4 
42.1 
42.1 
42.1 
42.  1 
42.1 


1  18 
1  19 
1  20 
1  21 


22 
23 


1  24 


42 
42 
42 
42 
42 
42.1 
3.8 


54.37 
53.35 
50.53 
47.58 
44.51 
41.35 
38.12 
34.83 
31.49 
28.11 
24.69 
21.25 
19.99 


263.10 
264.76 
268.75 
272.27 
275.37 
278.14 
280.62 
282.88 
284.94 
286.84 
288.61 
290.27 
290.85 


450 

449 

448 

447 

446 

445 

443 

442 

441 

440 

438.9 

437.9 

437.5 


0 

1 

5 

10 


00 
39 
55 
78 


17.70 
26.48 
32.24 
27.03 
18.10 
10.97 
5.61 
1.36 
-0.00 


340.48 

341.68 

345.74 

351.87 

2.15 

21.39 

55.33 

90.21 

110.45 

121.17 

127.48 

131.58 

132.74 


2436 

2285 

1892 

1514 

1167 

892 

773 

875 

1141 

1484 

1862 

2255 

2400 


69.23 
69.19 
68.57 
66.79 
63.08 
56.95 
52.35 
56.40 
62.69 
66.60 
68.51 
69.20 
69.25 


135.43 
136.92 
140.46 
143.50 
146.09 
148.32 
150.23 
151.86 
153.26 
154.46 
155.46 
156.31 
156.58 


Figure  11.  Output  from  ''times  above  horizon"  option.   At  =  1  minute.  Epoch  of  date. 


iterated  until  the  time  at  which  the  satellite  is  on  the  horizon.  All  of  these  times  and 
positions  are  then  output  to  the  screen.  Note  that  the  information  printed  out  in  this 
option,  with  the  exception  of  the  apogee/perigee  positions,  is  identical  to  that  printed 
out  in  the  "ephemeris"  option.  By  making  At  large,  say  10  minutes,  in  the  "times  above 
horizon"'  option,  then  only  the  times  of  satellite  rise,  culmination,  and  set  are  computed 
and  output  (see  Fig.  12). 

Each  line  in  the  output.dat  file  consists  of  the  station  identifier  and  the  satellite  iden- 
tifier followed  by  the  same  data,  and  in  the  same  order,  as  that  printed  to  the  screen.  This 
can  be  seen  by  comparing  the  first  six  lines  of  Fig.  12  with  the  same  lines  in  Fig.  9. 


Satellite:  11111 


Station:  001   Daisy,  Tenn. 


year  mo  da  hr  mn   sec    lat 


height  look 

long    (km)    elev    azim   range  angle   head 


1983  2  1 

1983  2  1 

1983  2  1 

1983  2  1 

1983  2  1 

1983  2  1 


1  13  19.4 
1  18  42.1 
1  24   3.8 


54.37  263.10 
38.12  280.62 
19.99  290.85 


2  49  40.9  45.38  250.95 
2  54  12.9  30.63  261.84 
2  58  44.7   15.05  269.44 


1983  2  1  14  47  13.6  16.76  285.33 
1983  2  1  14  50  58.0  29.64  291.62 
1983   2   1  14  54  43.3   42.02  300.04 


1983  2  1  16  20  48.3 

1983  2  1  16  26  16.7 

1983  2  1  16  31  46.6 

1983  2  1  17  59  48.3 

1983  2  1  18   2  51.4 

1983  2  1  18   5  55.6 


17.88  262.22 
36.50  272.26 

53.37  289.30 

37.38  249.28 
47.05  257.60 
55.73    269.75 


450.1  0.00    340.48  2436    69.23    135.43 

443.8  32.24      55.33  773    52.35    150.23 

437.5  -0.00    132.74  2400    69.25    156.58 

446.7  0.00    306.66  2428    69.23    145.41 
441.0  11.96   251.45  1429   66.15    153.58 

436.2  -0.00    195.12  2394    69.26    157.45 

433.6  0.00    150.29  2387   69.32      22.83 

437.9  6.44    106.44  1792   68.35      26.06 
443.2  -0.00      63.10  2419    69.28      32.12 

433.9  0.00   216.51  2389   69.32      23.02 

440.8  55.22   303.76  528   32.33      28.94 
448.2  -0.00      24.72  2431    69.26      43.11 

441.2  0.00    283.77  2414    69.29      29.38 

445.5  3.50    317.85  2066    69.01      36.02 

449.2  -0.00    351.87  2433    69.26      46.77 


Figure  12.  Output  from  "times  above  horizon"  option.   Ai  =  10  minutes.  Epoch  of  date. 

Fig.  13  is  a  sample  output  obtained  using  the  NAVSPASUR  reference  time.  A  compar- 
ison of  Fig.  10  and  Fig.  13  shows  that  the  satellite  latitude,  height,  and  heading  is  the 
same  but  the  longitude  differs  because  of  the  difference  of  the  earth's  rotation  between  the 
two  different  reference  times.  The  satellite-station  relationships  differ  for  the  same  rea- 
son. Figures  12  and  14  are  completely  dissimilar  because  of  the  differing  satellite-station 
relationships. 


Satellite:  11111 


Station:  001   Daisy,  Tenn. 


Apogee : 
Height 
Latitude 
Longitude 

Perigee : 
Height 
Latitude 
Longitude 


442.8   km. 

-55.1540 

230.8483 


450.4   km. 
55.1539 
50.8483 


Tue  1  Feb  1983 


year  mo  da  hr  mn   sec    lat 


height  look 

long     (km)    elev    azim   range  angle   head 


1983 
1983 
1983 
1983 
1983 
1983 
1983 


0   0   0.0  -12.26   98.20   434.0  -77.96  351.78  12919  11.09  157.82 


0  10  0.0  -45.84  118.32 

0  20  0.0  -65.19  179.03 

0  30  0.0  -44.71  237.62 

0  40  0.0  -10.95  257.18 

0  50  0.0   24.08  271.70 

1  0  0.0   55.55  299.61 


441.5  -79.40  232.53  12965  10.02  145.02 

446.4  -60.98  210.44  11620  27.15  88.46 

439.1  -41.32  205.69  9023  44.80  34.07 
430.7  -20.38  203.40  5458  61.42  22.04 
435.9   12.73  194.92   1373  65.77  24.38 

449.2  -5.07   32.07   3062  68.64  46.47 


Figure  13.  Output  from  ephemeris  option.  NAVSPASUR  reference  time. 


Satellite:    11111  Station:    001        Daisy,    Tenn  . 

height  look 

year   mo   da    hr   mil      sec         lat         long  (km)         elev        azim      range    angle      head 

1983  2  1  0  47  32.1  15.54  267.79  433.2  0.00  199.85  2385  69.32  22.63 
1983  2  1  0  52  59.8  34.24  277.36  439.8  58.70  112.56  509  29.01  27.87 
1983      2      1      0    58    28.6      51.41    293.03        447.4      -0.00      33.16      2429    69.27      40.56 

1983  2  1  2  24  57.1  29.86  251.15  438.0  0.00  262.00  2405  69.30  26.14 
1983  2  1  2  29  10.7  43.72  260.91  443.9  8.69  312.74  1641  67.65  33.31 
1983      2      1      2    33    25.7      55.99    276.83        449.3      -0.00        3.07      2433    69.26      47.22 

1983  2  1  8  58  59.6  56.01  272.80  450.7  0.00  356.74  2437  69.23  132.74 
1983  2  19  3  16.0  43.68  288.80  446.0  8.83  47.32  1637  67.55  146.72 
1983      2      19      7    31.7      29.71    298.62        440.7      -0.00      98.40      2412    69.24    153.92 

1983  2  1  10  33  57.2  51.41  256.65  449.0  0.00  326.74  2433  69.23  139.44 
1983  2  1  10  39  26.8  34.21  272.34  442.3  58.17  247.39  514  29.48  152.14 
1983      2      1    10   44    55.7      15.45    281.93        436.3      -0.00    160.33      2394    69.26    157.39 

1983  2  1  12  13  31.5  31.47  250.42  441.3  0.00  266.88  2414  69.24  153.27 
1983  2  1  12  14  27.8  28.30  252.20  440.2  0.29  256.89  2378  69.24  154.39 
1983      2      1    12    15    23.2      25.15    253.85        439.1      -0.00    247.03      2407    69.25    155.34 

Figure  14.  Output  from  "times  above  horizon"  option.  i\t  =  10  minutes.   NAVSPASUR  reference  tune. 

III.  SOURCES 

As  mentioned  in  the  introduction,  the  high  resolution  and  simplified  satellite  kinematics 
models  were  obtained  from  NAVSPASUR  [Ref.  1].  The  high  resolution  model  is  based  upon 
a  paper  by  Brouwer  [Ref.  2].  The  code  for  the  high  resolution  model  is  given  Appendix  E, 
while  the  code  for  the  simplified  model  is  contained  in  the  pos. update  subroutine  in 
Appendix  D.  Both  of  these  routines  require  a  solution  of  Kepler's  equation.  A  non-iteral 
state-of-the-art  solution  of  Kepler's  equation,  based  upon  a  paper  by  Mikkola  [Ref.  6]  is 
contained  in  subroutine  kepler  listed  in  Appendix  D. 

The  ground  track  determination  is  essentially  that  given  by  Escobal  [Ref.  3.  Ch.  10.3.2] 
but  with  a  minor  improvement  described  by  Shudde  [Ref.  7]. 

The  method  of  computation  of  the  times  of  satellite  culmination,  rise,  and  set  are  given 
in  App'  lix  B,  and  the  computation  of  the  satellite  heading  is  contained  in  Appendix  C. 
Additional  references  are  given  there. 

The  IAU  (1976)  System  of  Astronomical  Constants  [Ref.  10,  pg.  S7]  is  used  in  SATSTA 
as  well  as  the  J2000.0  reference  frame  [Ref.  10,  pg.  LI].  Satellite  positions  computed  by 
NAVSPASUR's  SHOWALL  [Ref.  1]  use  the  WGS  72  constants  and  the  J1900.0  reference  frame. 
For  this  reason,  the  agreement  between  SATSTA  and  SHOWALL  is  not  exact. 

IV.  DISCUSSION 

Although  the  position  updating  in  SATSTA  is  based  upon  the  satellite  kinematics  model 
of  Brouwer  [Ref.  2],  other  models  are  currently  in  use.   Most  notable  is  the  work  of  Hoots 


and  Roehrich  [Rcf.  5].  SATRAK  [Ref.  4]  is  a  FORTRAN  program  implementing  their  models 
and  is  available  from  the  Aerospace  Defense  Command,  Peterson  AFB.  Unfortunately,  the 
source  code  is  not  available  for  modification. 

The  program  presented  here  can  be  used  as  a  basis  for  a  more  sophisticated  model. 
The  input  scheme  could  be  improved  by  using  an  input  editor  to  catch  improper  types  of 
entry  and  to  allow  corrections  or  changes  to  be  made  to  previous  entries.  File  management 
programs  can  be  used  to  select  satellite  elements  and  station  parameters  from  master  files 
and  insert  them  into  "working"  files  for  program  input. 

Various  sensor  types  have  not  been  modelled  and  no  provision  has  been  made  for 
including  them  in  the  program.  Also,  satellite-satellite  visablity  has  not  been  modelled.  If 
such  features  are  desirable,  they  could  be  included  by  modifying  the  existing  program. 
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APPENDIX  A.  THE  SIMPLIFIED  KINEMATICS  ALGORITHM. 

The  computation  of  the  perturbed  satellite  orbit  about  the  earth  is  performed  using  a 
method  which  is  an  approximation  of  that  used  by  NAVSPASUR  [Ref.  1].  According  to 
NAVSPASUR,  this  approximation  will  yield  satellite  positons  with  a  random  error  of  about 
10-15  nautical  miles  for  a  20  day  interval  around  epoch. 

The  mean  elements  at  epoch,  which  must  be  input,  are: 

Eccentricity  eq 

Right  ascension  of  the  ascending  node  Q0 

Inclination  iQ 

Argument  of  perigee  u?o 

Mean  anomaly  Mq 

Mean  motion  Mi,  and 

Decay  M2 
Using  the  mean  elements,  the  following  auxiliary  quantities  are  computed: 

a0  =  m;2,\ 

r  1  ■       .  1 2/3 

C(3  cos    iq  —  1 
a0  =  a0    1  +  — — — 

5  cos2  z'o  —  1 
al(l-el) 

ao\l      eo) 
4a0 


(A-la) 


■•>  =  C    2n     "2^2,  (A-2) 


a  = 


3A/i 


M2,  (A-4) 


where  C  =  3/4^2  is  related  to  the  oblateness  of  the  earth.  NAVSPASUR  states  the  following: 

For  optimal  accuracy,  ao  from  Equation  A-l  is  used  in  Equation  A-la  to  obtain  a  final 
value  for  ao.  This  latter  ao  is  then  used  throughout.  Where  program  size  is  a  major 
consideration  and  additional  errors  of  some  ten  miles  in  satellite  height  can  be  tolerated 
the  user  may  choose  to  leave  out  Equation  A-la  and/or  Equation  A-4. 

Using  the  mean  values  of  the  elements  at  epoch  To  and  the  results  of  the  above 

calculations,  the  mean  values  of  the  elements  for  any  time  t  are: 

M  =  M0  +  Mi(t  -  To)  +  M2{t  -  T0)\ 
e  =  c0, 

u  =  ldq  +  Cj{M  -  Mq), 
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n  =  Qo  +  Cl(M  -  M0)  -  ue(t  -  r0), 
i  =  Zo ,      and 
a  =  a0  +  a(t  —  T0), 


where  u;e  is  the  earth's  sidereal  rotation  rate. 
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APPENDIX  B.  THE  TIME  OF  CULMINATION,  RISE  AND  SET  ALGO- 
RITHM. 

In  this  program  the  earth  is  considered  to  be  non-rotating,  i.e.,  the  coordinates  of  a  station 
are  given  by  latitude  and  longitude  with  the  longitude  unaffected  by  the  earth's  rotation. 
Rotation  is  accomplished  by  subtracting  the  accrued  rotation  from  the  right  ascension  of 
the  ascending  node  of  the  satellite.  Although  this  departs  from  realism,  the  benefit  is  that 
the  rectangular  coordinates  of  the  station  need  to  be  calculated  only  once. 


Fig.  B-l.  Satellite/Station  Geometry. 


Referring  to  Fig.  B-l,  culmination  of  the  satellite,  S,  with  respect  to  a  ground  station. 
P,  occurs  when  the  elevation  angle,  h,  of  the  satellite  is  a  maximum,  which  is  equivalent 
to  the  minimization  of  the  of  the  zenith  angle,  90°  —  h,  which  is  also  equivalent  to  the 
minimization  of  the  angle  77  between  z,  the  unit  vector  normal  to  the  station  location,  and 
r,  the  vector  from  the  earth's  center,  0,  to  the  satellite.  The  equation  for  i)  is 


z    r 


cos  77  = 


(B-l) 


where  z  is  given  by 


z  = 


cos  4>  cos  A 

cos  0sin  A 

sin</> 


(B-2) 


and  4>  and  A  are  the  station  geodetic  latitude  and  longitude,  respectively.    The  satellite 
position  vector  r  is  obtained  from  [Ref.  3,  Equ.  (3.57)] 


r  =  x^P  +  y^Q. 
13 


(B-3) 


where  xw  and  yw  measure  the  satellite  position  in  the  orbital  plane.1  P  and  Q  are  the 
first  two  columns  of  the  Euler  matrix  [Ref.  3,  Equs.  (3.43)  and  (3.44)] — they  transform  the 
in-plane  xu  and  yw  coordinates  to  the  equatorial  rectangular  system  and  are  given  by 


P  = 


and 


Q  = 


Qy 


cos  (jj  cos  r  —  sin  u  sin  T  cos  i 

cos  u)  sin  T  +  sin  u  cos  T  cos  i 

sin  u)  sin  i 

—  sin  uj  cos  T  —  cos  w  sin  T  cos  i 

—  sin  u;  sin  T  +  cos  w  cos  T  cos  z 

coso;  sinz 


(B-4a) 


(B-4b) 


where  u  is  the  argument  of  perigee  and  i  is  the  orbital  inclination;  T  =  Q,  —  u>e(t  —  tQ), 
where  Q  is  the  longitude  of  the  ascending  node  at  the  time,  i0,  of  the  last  element  update; 
ue  is  the  earth's  sidereal  rotation  rate  and  t  is  the  current  time. 

The  components  of  r  vary  with  time.  We  will  transform  the  equations  in  such  a 
way  that  will  replace  time  by  the  eccentric  anomaly  E  as  the  independent  variable.  The 
inclination  i  is  time  independent  and  the  elements  0,  and  u  will  be  considered  to  be 
independent  of  time  over  the  period  of  a  single  orbit  even  though  they  are  time  dependent 
through  perturbations.  The  components  of  r  which  will  be  considered  to  be  time  dependent 
are  then  xu,  yw  and  T. 

Write  Equation  (B-l)  as 

z  •  r  —  r  cost/  =  0. 


Take  the  derivative  with  respect  to  E  to  obtain, 

dr        dr 


dr\ 


dE 


_,  cos  77  +  r  sin  rj  —  =  0. 
dE  dE 


To  find  the  value  of  E  which  minimizes  77,  set  dr\jdE  =  0.  Thus, 


dr        dr 
ZdE-dECOSV  =  °' 


or 


dr        1  dr 


:Z-r  =  0. 


dE       r  dE 
Since  this  equation  cannot  be  solved  explicitly  for  E,  set 

dr        1  dr 


F(E) 


dE       r  dE 


z    r 


(B-5) 


The  origin  of  the  system  is  at  the  focus  of  the  orbit.  Jw  is  positive  in  the  direction  of  the  perifocus.  yu 
lies  in  the  orbital  plane  and  is  orthogonal  to  Xu. 
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and  find  E  such  that  F(E)  =0  by  the  Newton- Raphson  method.   The  Newton- Raphson 
procedure  also  requires  F'(E),  the  derivative  of  F(E)  with  respect  to  E: 

„.,„,  d2r         1    /  dr\2       1  d2r  1  dr         dr 

The  first  and  second  derivatives  of  Equation  (B-3)  are 

dr        dxu  dP       dy^  dQ 


and 


d2r 
dE2 

= 

d2xu          ndx^dP 
dE2             dE  dE        L 

d2P 
"dE2 

+ 

d2y^ 
dE2 

where 

[Ref. 

3. 

Equ.  (3.69)] 

dy^dQ  d2Q 


x^  =  a(cos  E  —  e)     and 

y«j  —  a  V  1  —  e2  sin  E. 

with  first  and  second  derivatives 

dx^  . 

=  —asm  is. 


=  a  v  1  —  e2  cos  E, 
an, 

and 


dE 

dy^ 

dE 

rf2y. 

i 

dE2 

d2x^ 

i 

a  v  1  -  e2  sin  £  =  — yw, 

rf£2    =-acos£. 
The  first  and  second  derivatives  of  Equations  (B-4)  axe 

dPT  _  oT 


r 


=  -^T^, 


(B-Sa) 


B-Sb) 


(B-Sc) 


dE  vdE 

dE         xdE' 
^  =  0, 

d£  ~Vrd£" 

— —  =  o 

dE         ' 
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an( 


d2px 

dE2 

=  -Px| 

(dY\ 

\dE) 

\-P*T 
1          ydE2' 

d2Py 

dE2 

=  -Py\ 

fdT\ 
KdE) 

2     p  rf2r 

d2Pz 
dE2 

=  0, 

d2Qx 
dE2 

=  -Qx 

\dE, 

I        WydE2' 

d2Qy 
dE2 

=  ~Qy 

\dEj 

)    +  Q*dE» 

d2Qz 

-  n 

(B-9b) 


and 


dE2 
Also  required  is  [Ref.  3,  Equ.  (3.67)] 

r  =  a(l  -ecosE)  (B-lOaJ 

and  it's  derivatives 


B-lOb) 


dr 

— —  =  aesini?,      and 
dE 

d2r 

—  =aecosE. 

Finally,  using  Kepler's  equation  [Ref.  3,  Equ.  (3.79)], 

E  —  e  sin  E      ^  _: 

*  = +T,  B-ll) 

n 

where  n  is  the  rate  of  orbital  motion  and  T  is  the  time  of  last  perifocal  passage,  the 

equation  for  T  can  be  written  as 

r  =  n-UJE-e™E+T-tX  (B-12a) 

Then, 

s-=^.  (B.12c) 

To  solve  Equation  (B-5),  an  initial  estimate  of  E  is  made  and  the  right-hand  side 
of  Equation  (B-13)  is  evaluated.  The  result,  the  left-hand  side  of  Equation  (B-13),  is  an 
improved  estimate  of  E, 

E  =  E-JWY  (B"13) 
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This  iteration  is  continued  until  \F(E)/F'(E)\  is  less  than  some  prescribed  amount.  The 
value  of  E  obtained  will  either  maximize  or  minimize  77.  The  desired  value  of  r\  will  be  the 
one  for  which  cos  77  >  0. 

Once  the  eccentric  anomaly  for  culmination,  Ec,  is  found,  the  time  of  culmination  is 
determined  directly  from  Equation  (B-ll).  Perturbing  Ec  by  plus  or  minus  about  10°,  it 
is  possible  to  obtain  good  starting  values  for  the  determination  of  the  satellite  set  and  rise 
times,  respectively. 

The  procedure  for  finding  the  rise  and  set  times  is  very  similar  to  finding  the  time 

of  culmination.    The  rise/set  equations  are  developed  by  Escobal  [Ref.  3,  Sec.  5.4],  but 

the  algorithm  presented  here  is  somewhat  simplified  since,  unlike  Escobal,  the  time  of 

culmination  is  available.    In  fact,  many  of  the  same  equations  can  be  used.    Again,  the 

eccentric  anomaly  will  be  used  as  the  independent  variable.  Referring  to  Fig.  B-l,  satellite 

rise  or  set  occurs  when  the  elevation  angle  h  is  equal  to  zero,  where  h  can  be  found  from 

the  equation 

smh  =  — ,  (B-14) 

P 

where 

p  =  r-R.  (B-15) 

In  Equation  (B-15),  R  is  the  station  coordinate  vector  and  p  is  the  slant  range  vector 
from  the  station  to  the  satellite.  For  the  spheroid  earth  [Ref.  3,  Equ.  (5.54)], 


R 


G\  cos  <f>  cos  A 

G\  cos  <f>  sin  A 

G2  sin  <p 


B-16' 


where  (f>  and  A  are  the  station  geodetic  latitude  and  longitude,  respectively.    G\  and  G2 
are  factors  which  account  for  the  earth's  oblateness,  and  are  given  by 


Gl  =  e  +  H         and 

V/l-(2/-/2)sin20 

G2=     ,      (1-/2)"  +* 


(B-17) 


v/l-(2/-/2)sin2^ 


where 


ae  =  earth's  equatorial  radius, 
/  =  earth's  flattening  factor,  and 
H  =  station  elevation  above  the  ellipsoid. 
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Write  Equation  (B-14)  as 

p  ■  z  =  psin  h. 

When  the  satellite  is  on  the  horizon,  h  =  0  and  sin  h  =  0.  Since  the  resulting  equation, 
p  •  z  =  0,  cannot  be  solved  explicitly  for  E,  set 

G(E)  =  p-z 

and  find  E  such  that  G(E)  =  0  by  the  Newton-Raphson  method.  Using  Equation  (B-15), 
write  G(E)  as 

G{E)  =  z  •  (r  -  R).  (B-18) 

The  Newton-Raphson  procedure  also  requires  G'(E),  the  derivative  of  G(E)  with  respect 
to  E1.  Since  the  earth  is  being  treated  as  non-rotating  with  the  sidereal  rotation  subtracted 
from  the  satellite's  ascending  node,  R  is  independent  of  E,  hence 

G'(E)  =  *~.  (B-19) 

In  Equations  (B-18)  and  (B-19),  z,  r,  dr/dE,  and  the  additional  needed  derivatives  have 
all  been  given  previously  in  this  section. 

To  solve  Equation  (B-18),  and  initial  estimate  of  E  is  made  and  the  right-hand  side 
of  Equation  (B-20)  is  evaluated.  The  result,  the  left-hand  side  of  Equation  (B-20),  is  an 
improved  estimate  of  £", 

E  =  E-§WY  (B"20) 

The  iteration  is  continued  until  \G(E)/G'(E)\  is  less  than  some  prescribed  amount.  The 
corresponding  time  of  rise  or  set  is  then  obtained  from  Equation  (B-ll). 
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APPENDIX  C.  THE  SATELLITE  HEADING  COMPUTATION. 

A  satellite  can  be  located  by  a  position  vector  x  =  T]i  +  x2j  +  x3k  in  a  rectangular 
coordinate  system  with  origin  at  the  earth's  center.  In  matrix  form,  and  with  Xj,  i2  and 
£3  expressed  in  spherical  coordinates, 


x  = 


r  cos  (f>  cos  A 

r  cos  <f>  sin  A 

r  sin  <f> 


where  <j>  is  the  latitude,  A  is  the  longitude  at  the  point,  and  r  =  |x|  is  the  distance  of  the 
point  from  the  earth's  center. 

Let  z  be  the  unit  vector  in  the  direction  of  r.  Then 


x 

Si 


-1 

z2 
Z3 


cos  4>  cos  A 
cos  <p  sin  A 

sin  o 


It  is  convenient  to  define  two  other  unit  vectors  in  the  plane  tangent  to  z.  These  are  local 
north  n  and  local  east  e,  defined  by 


n  = 


where 


—  sin  0cos  A 

—  sin  <£>sin  A 

cos  (f) 

dz 

z,  =  —  and 


e  = 


sin  A 
cos  A 

0 


za 


dz 


The  vectors  z,  n  and  e  provide  the  basis  for  a  right-handed  orthogonal  coordinate  system. 

If  a  satellite's  velocity  vector  v  =  v\\  +  v2]  +  t'3k  is  known  at  position  x  then  it  is 

easily  shown  [Ref.  8]  that 

n  •  v  =  cos  q     and 

e    v  =  sin  a, 

where  a  is  the  course  angle  or  heading.  From  these  relations,  a  is  found  to  be 

a  —  qatn(e  •  v,  n  •  v) 

where  qatn(y,x)  is  arctan(y/x)  corrected  to  the  proper  quadrant. 

Since  the  position  and  velocity  component's  are  in  rectangular  coordinates,  it  will  be 
necessary  to  convert  them  to  spherical  coordinates  in  order  to  determine  the  components 
of  n  and  e.  However,  the  spherical  coordinate  step  can  be  circumvented  as  follows:  First, 

cos  a  =  n  •  v  =  —  v\  sin  <p  cos  A  —  v2  sin  4>  sin  A  +  V3  cos  <t>. 
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Multiply  both  sides  by  cos  <f>  to  give 

cos  <f>  cos  a  =  —  i'i  sin  <j)  cos  <f>  cos  A  —  v2  sin  (j)  cos  <f>  sin  A  +  u3  cos2  0 

=  —  sin  0(2:1  vi  +  22^2)  +  V3  cos2  4>  +  Z3U3  sin<jf>  —  23U3  sin0 
=  -  sin  <(>(z\V\  +  z2v2  +  23V3)  +  i>3  cos2  0  +  23U3  sin  0 
=  -23(z  •  v)  +  t>3  cos2  0  +  v3  sin2  0 
=  ^3  -z3(z- v). 


Similarly, 


so  that 


Then, 


sin  a  =  e  •  v  =  — Vi  sin  A  +  v 2  cos  A, 
cos  <f>  sin  q  =  v2  cos  0  cos  A  —  V\  cos  0  sin  A 

=  Z\V2  —  2-2^1- 

a  =  qatn(cos  0  sin  a,  cos  0  cos  a). 


Since  cos  0  >  0  for  all  <p  £  (— 7r/2,7r/2),  the  cos  0  term  effectively  "cancels  out"  [Ref.  9] 
and  so  the  heading  is  determined  by 

a  =  qatn(sina,  cosq) 

=  qatn(zii>2  -  z2vx,v3  -  z3(z  •  v)). 

This  equation  is  used  by  NAVSPASUR  [Ref.  1]  in  SHOWALL.  In  SATSTA,  the  equivalent  form, 

(xxv2-x2vx             x3(x-v) 
a  =  qatn  ( ,  v3 


r' 


is  used. 


20 


APPENDIX  D.  THE  SATSTA  PROGRAM  LISTING. 

'SATSTA 

'  Satellite  Track.   03-03-88.         Rev.  09-26-89  0  0830. 

DECLARE  FUNCTION  acost  (xt) 

DECLARE  FUNCTION  asint  (xt) 

DECLARE  FUNCTION  atan2t  (yt,  xt) 

DECLARE  FUNCTION  decimalt  (v$) 

DECLARE  FUNCTION  dmodt  (xt,  mt) 

DECLARE  FUNCTION  dott  (at(),  bt()) 

DECLARE  FUNCTION  grastt  (tut,  tmt) 

DECLARE  SUB  apo. peri. gee  (elemt(),  elemut()) 

DECLARE  SUB  culmnate  (zvt(),  wet,  tOt,  eat,  cosetat,  eleraut()) 

DECLARE  SUB  euler  (perit,  nodet ,  inclt,  psnt(),  velt()) 

DECLARE  SUB  ground. track  (rt,  sat.dct,  sat.ltt,  sat.htt) 

DECLARE  SUB  hr.min.sec  (time  .of  .dayt ,  hr%,    minX,  sect,  nX) 

DECLARE  SUB  jd. to. date  (jdt(),  mk ,   dk ,    y*,  ht,  mo$,  day$) 

DECLARE  SUB  Julian . day .number  (m*.  dft ,  yft ,  utt,  jdt(),  djt,  tjt) 

DECLARE  SUB  kepler  (mt,  ect,    eat) 

DECLARE  SUB  output4  (k*/.,  yrft,  molt,  daft,  time .  of  .dayt ,  ltt,  lnt,  htt,  elt, 

azt,  rngt,  langlet,  headt,  sta.id$,  sat.id$) 
DECLARE  SUB  pos. update  (IX,  elemt(),  timet,  sat.xt(),  sat.xdt(),  elemu«(), 

rt ,  eat) 
DECLARE  SUB  rise. set  (zvt(),  sta.Rt(),  wet,  tOt,  eat,  elemutO) 
DECLARE  SUB  sat. head  (rt,  sat.xtQ,  sat.xdt(),  headt) 

DECLARE  SUB  sat .  sphere .  coord  (sat.xtQ,  rt,  tjt,  tmt,  sat.dct,  sat. lnt,  sat. rat) 
DECLARE  SUB  sta.p.coord  (ltt,  lnt,  htt,  xt()) 

DECLARE  SUB  sta.sat  (sta.xt(),  sat.xt(),  rngt,  azt,  elt,  langlet) 
DECLARE  SUB  time. update  (k'/.,  deltatt,  tmt,  jdt,  jdt(),  tjt,  time  .of  .dayt , 

mk ,  d& ,  yft,  ht ,  mo$,  day$) 
DECLARE  SUB  yrmoday  (v$,  yrft,  moft ,  dyft) 

'COLOR  7,  1 
'SCREEN  0,  1 

DEFDBL  A-Z 

OPTION  BASE  1 

DIM  jdepoch(2),  jdstart(2),  jdend(2) ,  jd(2)  ,  sat.x(3),  sat.xd(3),  sta.x(3,  4) 

DIM  elem(9),  elemu(9),  pv(3),  qv(3),  sta.z(3),  sta.R(3),  jd . charlie(2) 

'  constants 

CONST  pi  =  3.  141592653589793t,  rad  =  pi  /  180t,  twopi  =  pi  +  pi 

CONST  j2  =  .0010826318t 

CONST  nmi.ft  =  .3048t  /  1852t  '  n. mi/foot 

CONST  ae  =  6378. 14t  '  earth's  equatorial  radius  (km.) 

CONST  GE  =  398600. 5t  '  geo .  grav .  const,  (km) "3/(sec) "2 

k  =  60t  *  SQR(GE  /  ae)  /  ae'  (eru)~(3/2)/mm 

Herg.min  =  k 

CONST  fl  =  It  /  298.257t  '  earth's  flattening  factor 

CONST  e.sqr  =  (2t  -  fl)  *  fl  '  earth's  ellipticity  squared 

CONST  we  =  1.002737909350795t  *  twopi  /  1440t ' earth ' s  rotation  (radians/minute) 
CONST  sixty  =  60t 

'  Sample  input  (NORAD  Project  Spacetrack  Rpt .  t3 ,  Dec.  1980): 
amajor$  =  "1.04050189":  ecc$  =  "0.0086731":  aincl$  =  "72.8435" 
peri$  =  "52.6988":  anode$  =  "115.9689":  manomaly$  =  "110.5714" 
motion$  =  "16.05824518":  decay$  =  "0" 
epochdate$  =  "1980.1001":  epochtime$  =  "23.41241138" 
startdate$  =  "1980.1001":  atarttime$  =  "23.41241138" 
enddate$  =  "1980.1003":  endtime$  =  "0000" 
deltatS  =  "4" 

'  Sample  input  of  arbitrary  higher  altitude,  higher  ecc  orbit 

ecc$  =  "0.0200000":  aincl$  =  "64.94868" 

peri$  =  "202.277268":  anode$  =  "55.552464":  manomaly$  =  "126.670896" 
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motion!  =  "12.00000":  decay$  =  "4.281098d-3" 

epochdate*  =  "1983.0201":  epochtime$  =  "00.00" 

startdate$  =  "1983.0201":  starttime$  =  "00.00" 

enddate$  =  "1983.0208":  endtime$  =  "00.00" 

deltat$  =  "15" 

sta.lat$  =  "35.12":  sta.lon$  =  "-85.12":  sta.alt$  =  "500.0" 

'  Sample  input  (NAVSPASUR:  Sat.  22222  data  from  SHOWALL,  private  com. 

'     Mr.  Fred  Lipp,  July  1988): 

sat.id$  =  "22222" 

ecc$  =  "0.0023064":  aincl$  =  "64.94868" 

peril  =  "202.277268":  anode$  =  "55.552464":  manomaly$  =  "126.670896" 

motion$  =  "16.06302":  decay$  =  "4.281119d-3" 

epochdate$  =  "1983.0201":  epochtime$  =  "00.00" 

startdate$  =  "1983.0201":  starttime$  =  "00.00" 

enddate$  =  "1983.0203":  endtime$  =  "00.00" 

deltat$  =  "01" 

sta.id$  =  "001":  sta.name$  =  "Daisy,  Tenn." 

sta.latl  =  "35.12":  sta.lon$  =  "-85.12":  sta.altl  =  "500.0" 

'  Sample  input  (NAVSPASUR:  Sat.  11111  data  from  SHOWALL,  private  comm. 

'     Mr.  Fred  Lipp,  July  1988): 

sat.id*  =  "11111" 

ecc$  =  "0.0005545":  aincl$  =  "65.06057" 

peri$  =  "295.41470":  anode$  =  "272.43497":  manomaly$  =  "258.10682" 

motion!  =  "15.44194":  decay$  =  "0" 

epochdate$  =  "1983.0201":  epochtime$  =  "00.00" 

startdate$  =  "1983.0201":  starttime$  =  "00.00" 

enddateS  =  "1983.0202":  endtime$  =  "00.00" 

deltat$  =  "10" 

sta.id$  =  "001":  sta.name$  =  "Daisy,  Tenn." 

sta.lat$  =  "35.12":  sta.lon$  =  "-85.12":  sta.alt$  =  "500.0" 

start : 

CLS 

PRINT  "Select  data  input  mode  1,  2,  3  or  4:" 

PRINT 

PRINT  SPC(5);  "Keyboard  input  for  satellite  k   ground  station:" 

PRINT  SPC(IO);  "1.  Epoch  of  date." 

PRINT  SPC(IO);  "2.  NAVSPASUR  Reference  Date." 

PRINT 

PRINT  SPC(5);  "Data  file  input  for  satellite  &  ground  station:" 

PRINT  SPC(IO);  "3.  Epoch  of  date." 

PRINT  SPC(IO);  "4.  NAVSPASUR  One-Line  Charlie." 

data$  =  "" 

WHILE  data$  <  "1"  OR  data$  >  "4" 

data$  =  INPUT$(1) 
WEND 

CLS 

PRINT  "Select  run  mode  1  or  2:" 

PRINT 

PRINT  SPC(5);  "1.  Create  an  ephemeris." 

PRINT  SPC(5);  "2.  Find  times  satellite  is  above  the  horizon." 

case$  =  "" 

WHILE  case$  <>  "1"  AND  case$  <>  "2" 

case$  =  INPUT$(1) 
WEND 

GOSUB  input . sim. times 

SELECT  CASE  data$ 
CASE  "1",  "2" 

CLS 

PRINT  "Satellite  Data:" 

PRINT 
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PRINT  "Epoch  date  (yyyy.mradd)      ?  ";  epochdate$ 
LOCATE  3,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  epochdate* 
epochdate$  =  v$ :  CALL  yrmoday(v$,  yri,  mot,  dyfc) 

PRINT  "Epoch  time  (hh.mmss)        ?  ";  epochtime$ 

LOCATE  4,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  epochtime$ 

epochtime$  =  v$:  epochtime  =  decimalt(v$) 

CALL  Julian. day .number(mo*,  dy*,  yr*,  epochtime,  jdepoch(),  epoch,  tjepoch) 

LOCATE  6,  1:  PRINT  "Input  mean  elements  of  epoch:" 

PRINT  "Satellite  ID  Number         ?  ";  sat.id$ 
LOCATE  7,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  sat.id$ 
sat.id$  =  v$ 

PRINT  "Eccentricity  ?  ";  ecc$ 

LOCATE  8,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  ecc$ 
ecc$  =  v$:  ecc  =  VAL(v$) 

PRINT  "Long,  ascending  node  (deg)   ?  ";  anode$ 
LOCATE  9,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  anode$ 
anode$  =  v$:  anode  =  VAL(v$)  *  rad 

PRINT  "Inclination  (deg)  ?  ";  aincl$ 

LOCATE  10,  29:  INPUT  v$:  IF  v$  =  ""  THEN  v$  =  aincl$ 
aincl$  =  v$:  aincl  =  VAL(v$)  *  rad 

PRINT  "Arg.  of  perigee  (deg)        ?  ";  peri$ 
LOCATE  11,  29:  INPUT  v$:  IF  v$  =  ""  THEN  v$  =  peri$ 
peri$  =  v$ :  peri  =  VAL(v$)  *  rad 

PRINT  "Mean  anomaly  (deg)  ?  ";  manomaly$ 

LOCATE  12,  29:  INPUT  v$:  IF  v$  =  ""  THEN  v$  =  manomaly$ 
manomaly$  =  v$:  manomaly  =  VAL(v$)  *  rad 

PRINT  "Mean  motion  (rev/day)       ?  ";  motion! 

LOCATE  13,  29:  INPUT  v$:  IF  v$  =  ""  THEN  v$  =  motionS 

motion!  =  v$:  motion  =  VAL(v$)  *  twopi  /  1440t'  convert  to  radians/minute 

PRINT  "Decay  (rev/day"2)  ?  ";  decay$ 

LOCATE  14,  29:  INPUT  v$:  IF  v$  =  ""  THEN  v$  =  decay$ 
'  convert  to  radians/minute"2 : 
decay$  =  v$:  decay  =  VAL(v$)  *  twopi  /  (14401)  "  2 

LOCATE  16,  1:  PRINT  "Press  any  key  to  continue." 
FOR  i'/.  =  1  TO  10:  c$  =  INKEY$:  NEXT 
c$  =  INPUT$(1) 

CLS 

LOCATE  1,  1:  PRINT  "Station  Coordinates:" 

PRINT 

PRINT  "Station  Identification  Number  ?  ";  sta.id$ 

LOCATE  3,  42:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  sta.id$ 
sta.id$  =  v$ 

PRINT  "Station  Location  or  Name  ?  ";  sta.name$ 

LOCATE  4,  42:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  sta.name$ 
sta.namei  =  v$:  sta.lat  =  VAL(v$)  *  rad 

PRINT  "Station  latitude   (dd.mmss,  South  minus)  ?  ";  sta.lat$ 
LOCATE  5,  42:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  sta.lat$ 
sta.lat$  =  v$:  sta.lat  =  VAL(v$)  *  rad 

PRINT  "Station  longitude  (ddd.mmss,  West  minus)  ?  ";  sta.lon$ 
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LOCATE  6,  42:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  sta.lon$ 
sta.lon!  =  v$ :  sta.lon  =  VAL(v$)  *  rad 

PRINT  "Station  altitude   (feet)  ?  ";  sta.alt$ 

LOCATE  7,  42:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  8ta.alt$ 
sta.alt!  =  v$:  sta.alt  =  VAL(v$) 

LOCATE  9,  1:  PRINT  "Press  any  key  to  continue." 
FOR  it   =  1  TO  10:  c$  =  INKEY$ :  NEXT 
c$  =  INPUT$(1) 

OPEN  "output.dat"  FOR  APPEND  AS  #3 

GOSUB  compute 

CASE  "3" 

OPEN  "station.dat"  FOR  INPUT  AS  #1 
OPEN  "output.dat"  FOR  APPEND  AS  #3 

DO  UNTIL  EOF(l) 

INPUT  #1,  sta.id$,  sta.lat$,  sta.lon!,  sta.alt$,  sta.name! 
sta.lat  =  VAL(sta.lat$)  *  rad 
sta.lon  =  VAL(sta.lonS)  *  rad 
sta.alt  =  VAL(sta.alt$) 

OPEN  "elements.dat"  FOR  INPUT  AS  #2 

DO  UNTIL  E0F(2) 

INPUT  #2,  sat.id$,  epochdate!,  epochtime$,  ecc$,  anode$,  aincl$ ,  . 
peri$,  manomaly$,  motion$,  decay$ 

CALL  yrmoday(epochdate$,  yrJt ,  mot,  dyt) 

epochtime  =  decimal#(epochtime$) 

CALL  Julian .day .number(mo&,  dyt,  yrft,  epochtime,  jdepochO,  epoch, 

t jepoch) 
ecc  =  VAL(ecc$) 
anode  =  VAL(anode$)  *  rad 
aincl  =  VAL(aincl$)  *  rad 
peri  =  VAL(peri$)  *  rad 
manomaly  =  VAL(manomaly$)  *  rad 

motion  =  VAL(motion$)  *  twopi  /  1440#'  convert  to  radians/minute 
decay  =  VAL(decay$)  *  twopi  /  (1440#)  "  2   '  convert  to  rad/min"2 

jd(l)  =  jdstart(l) 
jd(2)  =  jdstart(2) 
jd  =  jd(l)  +  jd(2) 
time. of. day  =  starttime 

GOSUB  compute 

LOOP 
CLOSE  2 
LOOP 

CASE  "4" 

OPEN  "station.dat"  FOR  INPUT  AS  tl 
OPEN  "output.dat"  FOR  APPEND  AS  #3 

DO  UNTIL  E0F(1) 

INPUT  «1,  sta.id$,  sta.lat$,  sta.lon$,  sta.alt$,  sta.name$ 
sta.lat  =  VAL(sta.lat$)  *  rad 
sta.lon  =  VAL(sta. lon$)  *  rad 
sta.alt  =  VAL(sta.alt$) 

OPEN  "elements. nss"  FOR  INPUT  AS  «2 
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DO  UNTIL  E0F(2) 
INPUT  #2,  c$ 

sat.id*  =  MID$(c$,  1,  5) 

manomaly$  =  MID$(c$,  6,  9) 

motion$  =  HID$(c$,  14,  9) 

decay$  =  MID$(c$,  22,  7) 

ecc$  =  MID$(c$,  28,  9) 

pen$  =  MID$(c$,  36,  9) 

anode$  =  MID$(c$,  44,  9) 

aincl$  =  MID$(c$,  52,  9) 

epochdate$  =  "19"  ♦  MID$(c$,  60,  2)  ♦  "."  ♦  MID$(c$,  62,  4) 

epochtime$  =  "00.00" 

CALL  yrmoday(epochdate$,  yrft,  mot,  dyt) 

epochtime  =  decimal! (epochtime$) 

CALL  Julian. day. number (mot,  dy*,  yrt,  epochtime,  jdepoch(),  epoch, 

t jepoch) 
ecc  =  VAL(ecc$) 
anode  =  VAL(anode$)  *  twopi 
aincl  =  VAL(aincl$)  *  twopi 
peri  =  VAL(peri$)  *  twopi 
manomaly  =  VAL(manomaly$)  *  twopi 
motion  =  VAL(motion$)  *  Herg.min 
decay  =  VAL(.decay$) 

IF  decay  >  .51  THEN  decay  =  decay  -  If 
decay  =  decay  *  .00001  *  Herg.min  "  2 


convert  to  radians/minute 


convert  to  rad/min"2 


jd(l)  =  jdstart(l) 
jd(2)  =  jdstart(2) 
jd  =  jd(l)  +  jd(2) 

time. of. day  =  starttime 

G0SUB  compute 


LOOP 
CLOSE  2 

LOOP 
END  SELECT 
CLOSE 

END 


compute : 

'  initialize  times. 

tj  =  tjstart 

begin. time  =  (startdate  -  epoch)  *  1440f 

end. time  =  (enddate  -  epoch)  *  1440* 


'Julian  centuries  from  J2000.0 

'minutes 

'minutes 


elem(l 
elem(2 
elem(3 
elem(4 
elem(5 
elem( 6 
elem(7 
elem(8 
elem(9 


=  amajor  'semi-major  axis 

=  ecc  'eccentricity 

=  tzero  'time  of  perifocal  passage 

=  aincl  'inclination  (radians) 

=  anode  'longitude  of  ascending  node  (radians) 

=  peri  'argument  of  perifocus  (radians) 

=  motion  'mean  motion  (revolutions/day) 

=  manomaly  'mean  anomaly  (radians) 

=  decay  'decay  constant  (radians/minute"2) 


SELECT  CASE  data$ 
CASE  "1",  "3" 

'  Since  the  earth  is  treated  as  non-rotating,  elem(5)  must  be  corrected  for 

'   the  accrued  rotation  from  epoch  to  J2000.0 

elem(5)  =  elem(5)  -  15f  *  gmst(tj epoch,  0#)  *  rad 
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CASE  "2",  "4" 

'  NAVSPASUR  one-line  Charlie  correction 

CALL  Julian. day .number(l ,  1,  1985,  0,  jd. charlie() ,  dummy,  tj.charlie) 

vetc  =  (tjepoch  -  tj.charlie)  *  365251  -  (It  -  gmst(t j . charlie ,  Of)  /  24#) 
*  .9972695664t 

wetc  =  dmod(wetc  *  1440#  *  we,  twopi) 

elem(5)  =  elem(5)  +  wetc 

elem(5)  =  elem(5)  -  15t  *  gmst(tj epoch,  Ot)  *  rad 
END  SELECT 

CALL  sta.p.coord(sta.lat ,  sta.lon,  sta.alt,  sta.xQ) 
FOR  it   =  1  TO  3 

sta.z(iy.)  =  sta.x(i%,  2) 

sta.R(iy.)  =  sta.xCiy.,  1) 
NEXT 

CALL  pos .updated ,  elem(),  begin. time,  sat.x(),  sat.xd(),  elemu(),  r,  ea) 
CALL  pos  .update(2 ,  elem(),  begin. time,  sat.xQ,  sat.xdO,  elemu(),  r,  ea) 

CLS 

PRINT  "Satellite:  ";  sat.id$;  "     Station:  ";  sta.id$;  "    ";  sta.name$ 

PRINT 

SELECT  CASE  case$ 

> 

CASE  "1" 

tm  =  starttime  'UT  in  hours 

CALL  apo  .peri  .gee(elem() ,  elemuO) 

CALL  jd. to .date( jd() ,  mft,  dft,  y&,  ht,  mo$,  day$) 

PRINT 

PRINT  day$;  dft;  mo$ ;  yk 

CALL  output4(0,  yr& ,  mo&,  daft,  time. of. day,  sat.lt,  sat. In,  sat.ht,  el,  az , 
rng,  langle,  head,  sta.id$,  sat.id$) 

'  time-loop 

FOR  time  =  begin. time  TO  end. time  STEP  delta. t 

CALL  pos .  update(2 ,  elem(),  time,  sat.xQ,  sat.xdO,  elemu(),  r,  ea) 

'compute  spherical  coordinates 

CALL  sat . sphere . coord(sat .x() ,  r,  t j ,  tm,  sat.dc,  sat. In,  sat.ra) 

'satellite  heading 

CALL  sat.head(r,  sat.xQ,  sat.xdO,  head) 

'earth  ground  track 

CALL  ground. track(r,  sat.dc,  sat.lt,  sat.ht) 

'compute  station-satellite  relations 

CALL  sta.  sat(sta.x() ,  sat.xQ,  rng,  az,  el,  langle) 

rng  =  rng  *  ae 

CALL  output4(l,  yft,  mft,  dft,  time. of. day,  sat.lt,  sat. In,  sat.ht,  el,  az,  rng,  . 
langle,  head,  sta.idS,  sat.id$) 

'update  time 

CALL  time .update(0,  delta. t.hr,  tm,  jd,  jd(),  t j  ,  time. of. day,  mi,  dfc ,  y&,  _ 
hf,  mo$,  day$) 

NEXT  time 
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PRINT  :  PRINT  "Press  any  key  to  continue." 
FOR  i%   =  1  TO  10:  c$  =  INKEY$:  NEXT 
c$  =  INPUT$(1) 


CASE  "2" 

CALL  output4(0,  yt ,  m* ,  dft,  time . of .day ,  sat.lt,  sat. In,  sat.ht,  el,  az ,  rng , 
langle ,  head,  sta.id$,  sat.id$) 

tm  =  0 

period  =  twopi  /  motion 

ea  =  0 

time  =  begin. time 

initialize : 

CALL  pos .update(2 ,  elem(),  time,  sat.x(),  sat.xd(),  elemu(),  r,  ea) 

CALL  culmnate(sta.z() ,  we,  time,  ea,  coseta,  elemuO) 

IF  coseta  <  0  THEN 

ea  =  ea  ♦  pi 

CALL  culmnate(sta.z()  ,  we,  time,  ea,  coseta,  elemuO) 
END  IF 

time  =  (ea  -  elem(2)  *  SIN(ea))  /  motion  +  elemu(3) 
dt.hr  =  (time  -  begin. time)  /  sixty 
CALL  time .update(0 ,  dt.hr,  tm,  jd,  jd(),  tj ,  time. of. day,  m& ,  d& ,  y& ,  h# ,  _ 

mo$,  day$) 
IF  time  >  begin. time  ♦  period  THEN 

time  =  time  -  period 

CALL  time  .update(  1 ,  -period,  tm,  jd,  jd(),  t  j  ,  time. of. day,  m&  ,  die ,  y&,  . 
hf,  mo$,  day$) 

GOTO  initialize 
END  IF 

again : 

FOR  1%  =  1  TO  2 

CALL  pos  .update(2 ,  elem(),  time,  sat.xO,  sat.xdO,  elemuO,  r,  ea) 

CALL  culmnate(sta.z()  ,  we,  time,  ea,  coseta,  elemuO) 

dt.min  =  (ea  -  elem(2)  *  SIN(ea))  /  motion  ♦  elemu(3)  -  time 

time  =  time  +  dt.min 

CALL  time .updated ,  dt.min,  tm,  jd,  jd(),  t j ,  time. of. day,  mi,  dfe,  yft,  _ 
h«,  mo$ ,  day$) 
NEXT 
time.cul  =  time 

CALL  sta.sat(sta.x()  ,  sat.xO,  rng,  az,  el,  langle) 
IF  el  >  0  AND  time  >=  begin. time  THEN 

CALL  pos.update(2,  elem(),  time,  sat.xO,  sat.xdO,  elemuO,  r,  ea) 
ea  =  ea  -  .335 
FOR  i%  =  1  TO  2 

CALL  rise.set(sta.z() ,  sta.R(),  we,  time,  ea,  elemuO) 

dt.min  =  (ea  -  elem(2)  *  SIN(ea))  /  motion  ♦  elemu(3)  -  time 

time  =  time  +  dt.min 

CALL  time. updated ,  dt.min,  tm,  jd,  jd() ,  t j ,  time. of. day,  mft ,  dft,  _ 

yft,  hf ,  mo$,  day$) 
CALL  pos.update(2,  elem()  ,  time,  sat.xO,  sat.xdO,  elemuO,  r,  ea) 
NEXT 
time. rise  =  time 

CALL  8ta.sat(sta.xO,  sat.xO,  rng,  az,  el,  langle) 
CALL  sat . sphere .coord(sat.x() ,  r,  t j ,  tm,  sat.de,  sat. In,  sat.ra) 
CALL  ground. track(r,  sat.de,  sat.lt,  sat.ht) 
CALL  sat.head(r,  sat.xO,  sat.xdO,  head) 
CALL  jd.to.date(jd() ,  mft ,  dk,  y*.  ht,  mo$,  day$) 

CALL  output4(l,  yft,  m* ,  dft,  time. of. day,  sat.lt,  sat. In,  sat.ht,  el,  _ 
az ,  rng  *  ae ,  langle,  head,  Bta.idS,  sat.id$) 
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time.incr  =  dmod(time .cul  -  time. rise,  delta. t) 

time  =  time  ♦  time.incr 

CALL  time .updated ,  time.incr,  tm,  jd,  jd()  ,  t j ,  time. of. day,  mft,  dft, 

yft,  hi,  mo$,  day$) 
CALL  pos  .update(2,  elem(),  time,  sat.xO,  sat.xdO,  elemu(),  r,  ea) 
CALL  sta.  sat(sta.z() ,  sat.xO,  rng,  az,  el,  langle) 
CALL  sat . sphere . coord(sat .x() ,  r,  t j ,  tm,  sat.dc,  sat. In,  sat.ra) 
CALL  ground .track(r,  sat.dc,  sat.lt,  sat.ht) 
CALL  sat.head(r,  sat.xO,  sat.xdO,  head) 
CALL  jd.to.date(jd(),  mft,  dft,  yJt,  hi,  mo$,  day$) 

CALL  output4(l,  yft,  mft,  dft,  time. of. day,  sat.lt,  sat. In,  sat.ht,  el,  _ 
az ,  rng  *  ae ,  langle,  head,  sta.idl,  sat.id$) 

DO  WHILE  el  >  0 

time  =  time  +  delta. t 

CALL  time  .  update' 0 ,  delta. t.hr,  tm,  jd,  jd(),  t j ,  time. of. day,  mft,  dft, 
yft,  hi,  mo$,  day$) 

CALL  pos  .update(2 ,  elem(),  time,  sat.xO,  sat.xdO,  elemu(),  r,  ea) 

CALL  sta  .  sat (sta . x() ,  sat.xO,  rng,  az,  el,  langle) 

IF  el  <  0  THEN 
EXIT  DO 

END  IF 

CALL  sat .sphere . coord(sat . x() ,  r,  t j ,  tm,  sat.dc,  sat. In,  sat.ra) 

CALL  ground .track(r,  sat.dc,  sat.lt,  sat.ht) 

CALL  sat.head(r,  sat.xO,  sat.xdO,  head) 

CALL  jd. to.date( jd() ,  mft,  dft,  yft,  hf,  mo$ ,  day$) 

CALL  output4(l,  yft,  mft,  dft,  time. of. day,  sat.lt,  sat. In,  sat.ht,  el,  . 
az ,  rng  *  ae,  langle,  head,  sta.id$,  sat.id$) 
LOOP 

save. time  =  time 
FOR  i%  =  1  TO  2 

CALL  pos .  update(2 ,  elem() ,  time,  sat.xO,  sat.xdO,  elemu(),  r,  ea) 

CALL  rise .  set(sta.z()  ,  sta.RO,  we,  time,  ea,  elemuO) 

dt.min  =  (ea  -  elem(2)  *  SIN(ea))  /  motion  ♦  elemu(3)  -  time 

time  =  time  +  dt.min 

CALL  time .updated ,  dt.min,  tm,  jd,  jd(),  t j  ,  time. of. day,  mft,  dft,  yft, 
ht,  mo$,  day$) 
NEXT 

CALL  pos  .update(2 ,  elem(),  time,  sat.xO,  sat.xdO,  elemuO ,  r,  ea) 

CALL  sta.  sat(sta.x()  ,  sat.xO,  rng,  az ,  el,  langle) 

CALL  sat . sphere . coord(sat . x() ,  r,  t j ,  tm,  sat.dc,  sat. In,  sat.ra) 

CALL  ground .track(r,  sat.dc,  sat.lt,  sat.ht) 

CALL  sat.head(r,  sat.xO,  sat.xdO,  head) 

CALL  jd.to.date(jd() ,  mft,  dft,  yft,  hf,  mo$,  day$) 

CALL  output4(l,  yft,  mft,  dft,  time. of. day,  sat.lt,  sat. In,  sat.ht,  el,  _ 
az ,  rng  *  ae,  langle,  head,  sta.id$,  sat.idS) 

PRINT 

END  IF 

time  =  time. cul  ♦  period 
dt.hr  =  (time  -  begin. time)  /  sixty  -  tm 

CALL  time.update(0,  dt.hr,  tm,  jd,  jd(),  tj ,  time. of. day,  mft,  dft,  yft,  hf ,  _ 
mo$,  day$) 

IF  time  <  end. time  THEN  GOTO  again 
)__________________________________________________________________ 

END  SELECT 
RETURN 

input . sim. times : 
CLS 

PRINT  "Simulation-Time  Data:" 
PRINT 

PRINT  "Starting  date  (yyyy.mmdd)    ?  ";  startdate$ 
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LOCATE  3,  29:  INPUT  v$  :  IF  v$  =  ""  THEN  v$  =  startdateS 
startdate$  =  v$:  CALL  yrmoday(v$,  yr* ,  mo* ,  dyJt) 

PRINT  "Starting  time  (hh.mmss)      ?  ";  starttime$ 

LOCATE  4,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  starttime$ 

starttime!  =  v$ :  starttime  =  decimalf(v$) 

CALL  Julian,  day  .  number  ( mok  ,  dy*.  yrk,  starttime,  jdstartO,  startdate, 

t jstart) 
jd(l)  =  jdstart(l) 
jd(2)  =  jdstart(2) 
jd  =  jd(l)  ♦  jd(2) 
time. of. day  =  starttime 

PRINT  "Ending  date  (yyyy.mmdd)     ?  ";  enddate$ 
LOCATE  5,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  enddate$ 
enddate$  =  v$:  CALL  yrmoday(v$,  yr*.  mot,  dy*) 

PRINT  "Ending  time  (hh.mmss)       ?  "j  endtime$ 

LOCATE  6,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  endtime$ 

endtime$  =  v$ :  endtime  =  decimalf(v$) 

CALL  Julian .day . number(mot,  dyft ,  yrft,  endtime,  jdendQ,  enddate,  tjend) 

PRINT  "Time  increment  (minutes)    ?  ";  deltat$ 

LOCATE  7,  29:  INPUT  v$ :  IF  v$  =  ""  THEN  v$  =  deltatl 

deltat$  =  v$:  delta. t  =  VAL(v$) :  delta. t.hr  =  delta. t  /  sixty 

LOCATE  9,  1:  PRINT  "Press  any  key  to  continue." 
FOR  if.    =  1  TO  10:  c$  =  INKEY$ :  NEXT 
c$  =  INPUT$(1) 

RETURN 

FUNCTION  acosf  (x#)     '  02-10-88    Rev.  02-24-88 

'  The  arccosine  function  derived  from  the  ATN  function  with  no  singularities. 

'  The  angle  is  returned  in  the  (0,pi)  interval. 

'Note  that  in  relational  comparisons,  -1  is  "true"  and  0  is  "false". 

CONST  pi  =  3. 141592653589793t,  eps  =  1D-33 

IF  ABS(x#)  <=  If  THEN 

acost  =  ATN(SQR(1*  -  x#  *  x«)  /  (x«  -  eps  *  (xf  =  Of)))  -  pi  *  (xf  <  Of) 
ELSE 

PRINT 

PRINT  "Error  in  acos  function.  ABS(arg)  >  1" 

PRINT  "arg  =  ";  xf 

PRINT 
END  IF 
END  FUNCTION 

>**********^******* ***************************************************** *»»**** 

SUB  apo. peri. gee  (elemf(),  elemufQ)  '   09-14-88.     Rev  09-14-88. 

'  Compute  location  and  height  of  apogee  and  perigee  above  oblate  earth. 

'  NOTE:  (1)  The  output  from  this  routine  is  approximately  valid  for  the 

'  current  orbit  only.  In  particular,  the  longitude  is  not  corrected 

'  for  the  earth's  rotation  during  the  current  orbit.  Also,  the 

'  ascending  node  and  argument  of  perigee  are  not  corrected  for 

'  their  rates  of  change  during  the  current  orbit. 

'  (2)  For  nearly  circular  orbits,  the  height  of  the  satellite  will 

'  probably  get  loser  than  the  perigee  value  and  higher  than  the 

'  apogee  value  at  locations  other  than  the  perigee  and  apogee 

'  positions this  is  an  effect  of  the  earth's  oblateness. 
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CONST  pi  =  3.141592653589793t,  rad  =  pi  /  180t,  twopi  s  pi  +  pi 

'  spherical  earth  values: 

apogee!  =  elemut(l)  *  (It  -  elemt(2)) 

perigeei  =  elemut(l)  *  (It  ♦  elemt(2)) 

'  longitudes  and  declinations: 

sinclf  =  SIN(elem#(4)) 

spent  =  SIN(elemut(6)) 

apo.dct  =  asint(sperit  *  sinclf) 

apo.lnt  =  elemut(5)  ♦  atan2t(C0S(elemt(4))  *  sperit,  COS(elemuf (6) )) 

apo.lnf  =  dnodf (apo .lnf ,  twopi) 

IF  apo.lnt  >  pi  THEN  apo.lnt  =  apo.lnt  -  twopi 

'  next  are  approximations  should  update  by  peri-  ft  nodal-rates 

'     for  1/2  revolution: 
peri.dct  =  -apo.dct 
peri.lnt  =  dmod(apo.lnt  ♦  pi,  twopi) 
IF  peri.lnt  >  pi  THEN  peri.lnt  =  peri.lnt  -  twopi 

CALL  ground .track(apogeef ,  apo.dct,  apo.ltt,  apo.htt) 
CALL  ground. track(perigeef,  peri.dct,  peri.ltt,  peri.htt) 


PRINT 

PRINT  "Apogee:" 

PRINT  USING  "ft  tfft.t  ft" 

PRINT   USING    "ft   ttt.tttf 

PRINT  USING  "ft  !»f.ff*f" 

PRINT 

PRINT  "Perigee:" 

PRINT  USING  "ft  «*f«.f  ft" 

PRINT  USING  "ft  ttt.tttt" 

PRINT  USING  "ft  f#t.t#tt" 

PRINT 


"Height    =  ";  apo.htt;  "  km.  " 

"Latitude   =   ";  apo.ltt  /  rad 

"Longitude  =   ";  dmod(apo.lnt  /  rad,  360#) 


"Height     =  ";  peri.htt;  "  km.  " 

"Latitude   =   ";  peri.ltt  /  rad 

"Longitude  =   ";  dmod(peri .lnt  /  rad,  360f) 


IF  apo.htt  <  0  THEN 

PRINT  "WARNING:  Satellite  is  suborbital  and  will  impact  earth." 

PRINT  "Do  you  want  to  continue?  Yes  or  No." 

PRINT 

WHILE  c$  <>  "y"  AND  c$  <>  "n" 
c$  =  LCASE$(INPUT$(1)) 

WEND 

IF  c$  =  "n"  THEN  END 
END  IF 
END  SUB 

FUNCTION  asint  (x#)     '   02-10-88    Rev.  02-24-88 

'  The  arcsine  function  derived  from  the  ATN  function  with  no  singularities. 

'  The  angle  is  returned  in  the  (-pi, pi)  interval. 

'Note  that  in  relational  comparisons,  -1  is  "true"  and  0  is  "false". 

CONST  eps  =  1D-33 

IF  ABS(xt)  <=  It  THEN 

asint  =  ATN(xt  /  (SQR(lt  -  xt  *  xt)  -  eps  *  (ABS(xt)  =  It))) 
ELSE 

PRINT 

PRINT  "Error  in  asin  function.  ABS(arg)  >  1" 

PRINT  "arg  =  ";  xt 

PRINT 
END  IF 
END  FUNCTION 

>********»**»«*»*»**»»*******#*************»***»**♦*»**********«»»»****»**«»*** 
FUNCTION  atan2t  (yt,  xt)     '    02-10-88   Rev.  02-24-88 
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'  The  two  argument  quadrant  determining  arctangent  function. 
'  The  angle  iB  returned  in  the  (-pi, pi)  interval. 

'Note  that  in  relational  comparisons,  -1  is  "true"  and  0  is  "false". 

CONST  pi  =  3.141592653589793t,  epB  =  1D-33 

atan2t  =  ATN(yt  /  (it  -  eps  *  (xt  =  Of)))  -  pi  *  (xt  <  Of)  *  (SGN(yt)  _ 
-  (yt  =  Ot)) 

END  FUNCTION 

'ft***************************************************************************** 
SUB  culmnate  (zvf()  ,  vet,  tOt,  eat,  cosetaf,  elemut()) 

Compute  eccentric  anomaly  for  the  time  of  culmination  of  a  satellite 
01-30-89.  Rev.  02-02-89  C  1000. 

NOTE:  This  particular  form  of  the  culminate  function  was  written  from  the 

development  of  30  Jan  89.  It  assumes  a  non-rotating  earth  so  that  the 
station  longitude  is  fixed.  Rotation  is  accomplished  by  subtracting  the 
accumulated  rotation  from  the  right  ascension  of  the  ascending  node  of 
the  satellite.  This  is  the  method  used  to  compute  the  ephemeris  and  is 
also  equivalent  to  the  method  used  by  NAVSPASUR. 

lambda  =  station  geodetic  longitude 

phi  =  station  geodetic  latitude 

we  =  earth's  rotation  (sidereal  rate  of  change) 

tO  =  time  the  elements,  elemuQ,  were  last  updated 

capt  =  T  --  time  of  latest  perifocal  passage 

mot  =  mean  motion 

ea  =  eccentric  anomaly  E 

ecc  =  eccentricity  e 

amjr  =  semimajor  axis  (cancels  out  in  all  formulas,  so  not  included) 

CONST  pi  =  3. 141592653589793t,  rad  =  pi  /  180#,  twopi  =  pi  ♦  pi 

eccf  =  elemuf(2) 
captf  =  elemu#(3) 
inclf  =  elemui(4) 
nodef  -  elemuf(5) 
perif  =  elemuf(6) 
mott  =  elemuf(7) 

DIM  pv«(3),  pvdf(3),  pvddf(3),  qvt(3) ,  qvdt(3),  qvddf(3) 

DIM  rvect(3),  rvecdf (3) ,  rvecddf(3) 

zv    =  Z-vector ,  unit  normal  at  the  station 

pv    =  P-vector,  qv  =  Q-vector  (first  two  columns  of  the  Euler  matrix) 

rvec  =  R-vector  (satellite  position),  rvecd  =  dR/dE,  rvecdd  =  d2R/dE2 

rv    =  r magnitude  of  the  R-vector,  rvd  =  dr/dE,  rvdd  =  d2r/dE2 

gamma  =  node  -  accrued  rotational  motion 

coBwi  =  COS(perit) 
inwf  =  SIN(perif) 


s 

CO 


inw»  =  siKvperi»; 
osif  =  COS(inclf) 
inif  =  SIN(inclf) 


corrf  =  1 

WHILE  ABS(corrt)  >  .0001 

coseat  =  COS(eat) 

ecoseaf  =  eccf  *  coseat 

sineat  =  SIN(eat) 

esineaf  =  ecct  *  sineaf 

gammat  =  nodef  -  wet  *  ((eat  -  esineaf)  /  motf  ♦  captf  -  tOt) 
gammadf  =  -wet  *  (if  -  ecoseaf)  /  mott 
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gammaddt  =  -wet  *  esineaf  /  motf 

cosgf  =  COS(gammat) 
singf  =  SIN(gammaf) 

rtecct  =  SQR(lt  -  ecct  *  ecct) 

rv#  =  If  -  ecoBeat 
rvdt  =  esineaf 
rvddt  =  ecoseat 

zvt  =  coseaf  -  eccf 
xwdf  =  -sineat 
xwddt  =  -coseat 
ywt  =  rtecct  *  sineat 
ywdt  =  rtecct  *  coseat 
ywddt  =  -ywt 

pvt(l)  =  coswt  *  cosgt  -  sinwt  *  singt  *  cosit 
pvt(2)  =  coswt  *  singt  +  sinwt  *  cosgt  *  cosit 
pvt(3)  =  sinwt  *  sinit 

qvt(l)  =  -sinwt  *  cosgt  -  coswt  *  singt  *  cosit 
qvf(2)  =  -sinwt  *  singt  +  coswt  *  cosgt  *  cosit 
qvt(3)  =  coswt  *  sinit 

pvdt(l)  =  -pvt(2)  *  gammadt 
pvdf(2)  =  pvf(l)  *  gammadt 
pvdt(3)  =  Ot 

qvdt(l)  =  -qvf(2)  *  gammadt 
qvdt(2)  =  qvf(l)  *  gammadt 
pvdt(3)  =  Ot 

pvddf(l)  =  -pvdt(2)  *  gammadt  -  pvt(2)  *  gammaddt 
pvddt(2)  =  pvdf(l)  *  gammadt  ♦  pvt(2)  *  gammaddt 
pvddt(3)  =  Ot 

qvddt(l)  =  -qvdf(2)  *  gammadt  -  qvt(2)  *  gammaddt 
qvddf(2)  =  qvdf(l)  *  gammadt  +  qvt(2)  *  gammaddt 
pvddt(3)  =  Ot 

FOR  i%   =  1  TO  3 

rvect(i'/.)  =  xwt  *  pvt(i7.)  +  ywt  *  qvt(iX) 


rvectCi'/J  =  xwt  *  pvt(i7.)  +  ywt  *  qvt(i7.) 

rvecdt(iy.)  =  xwdt  *  pvt(iX)  ♦  xwt  *  pvdt(iy.)  +  ywdt  *  qvt(i'/.)  _ 

+  ywt  *  qvdt(i7.) 
tempt  =  xwddt  *  pvt(i7.)  ♦  2t  *  xwdt  *  pvdt(iy.)  +  xwt  *  pvddt(i'/.) 
rvecddt(i'/.)  =  tempt  +  ywddt  *  qvt(iX)  ♦  2t  *  ywdt  *  qvdt(i7.)  _ 
♦  ywt  *  qvddtUX) 
NEXT 

cosetat  =  dott(zvt(),  rvectQ)  /  rvt 

fet  =  dott(zvt(),  rvecdt())  -  rvdt  *  cosetat 

rrdt  =  rvdt  /  rvt 

fedt  =  dott(zvt(),  rvecddtO)  ♦  cosetat  *  (rrdt  -  rvddt) 

fedt  =  fedt  -  rrdt  *  dott(zvt(),  rvecdtO) 

corrt  =  fet  /  fedt 
IF  corrt  >  pi  THEN 

eat  =  eat  ♦  twopi 
ELSE 

eat  =  eat  -  corrt 

IF  cosetat  <  0  THEN  eat  =  eat  +  pi 
END  IF 


WEND 
END  SUB 
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FUNCTION  decimal!  (v$)  '   03-21-88.     Rev  03-21-88. 

'  Convert  ddd .mrassf f  or  hh.mmssfff  to  decimal  degrees  or  decimal  hours. 

ix  =  INSTR(v$,  ".") 
IF  ix  =  0  THEN 

decimal*  =  VAL(v$) 
ELSE 

xf  =  VAL(LEFT$(v$,  ix)) 

sn  =  0 

IF  xt  <  Ot  THEN 
sn  =  -1 
xt  =  -xt 

END  IF 

w$  =  v$  +  "0000" 

yt  =  VAL(MID$(w$,  ix  +  1,  2)) 

zt  =  VAL(MID$(w$,  ix  ♦  3,  2)  ♦  "."  ♦  RIGHT$(w$,  LEN(w$)  -  ix  -  4)) 

xt  =  (zt  /  60t  ♦  yt)  /  60t  ♦  xt 

IF  sn  THEN  xt  =  -xt 

decimalt  =  xt 
END  IF 
END  FUNCTION 

'*»»**«*»»*********»***»**»*»»******»*«***«**»»»*. ■...•••»............. ........ 

FUNCTION  dmodt  (xt,  mt)   '    02-10-88    Rev.  02-24-88 
'  Provides   x  MOD  m   for  real-valued  functions. 

IF  m«  <>  Ot  THEN 

dmodt  =  xt  -  mt  *  INT(xt  /  mt) 
ELSE 

PRINT 

PRINT  "Error  in  dmod  function.  Modulus  cannot  be  zero." 

PRINT 
END  IF 
END  FUNCTION 

FUNCTION  dott  (at(),  bt()) 
'  Dot  product  of  two  vectors. 

dott  =  at(l)  *  bt(l)  ♦  at(2)  *  b#(2)  ♦  at(3)  *  bt(3) 

END  FUNCTION 

<♦***«***.+*».**«***•***»*****»******»*»***»****»»**»»»**»«*»»»**»««*»».•»»»••» 
SUB  euler  (perit,  nodet,  inclt,  psnt(),  velt()) 
'  09-05-88.  Rev.  01-16-89. 

'Compute  the  first  two  columns  (p  k   q)  of  the  Euler  matrix. 

'Rotate  in-plane  position  and  velocity  components  to  rectangular  equatorial 

'    coordinates. 

'psn()  and  vel()  are  both  input  and  output  vectors.  It  is  assumed  that  the 

'    third  component  of  both  psn()  and  vel()  is  zero  on  input. 

cwt  =  COS(perit) 

Bwt  =  SIN(perit) 

cnt  =  COS(nodet) 

snt  =  SIN(nodet) 

cit  =  COS(inclt) 

sit  =  SIN(inclt) 

pxt  =  cwt  *  cnt  -  sw#  *  snt  *  cit 

pyt  =  cwt  *  snt  ♦  swt  *  cnl  *  cit 

pzt  =  swt  *  sit 

qxt  =  -swt  *  cnt  -  cwt  *  snt  *  cit 

qyt  =  -swt  *  snt  ♦  cwt  *  cnt  *  cit 

qzt  =  cwt  *  sit 
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* 
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* 

t2t 

* 

tit 

♦ 

qyt 

* 

t2t 

* 

tit 

+ 

qzt 

* 

t2t 

tit  =  psnt(l) 
t2t  =  psnt(2) 
psnt(l)  =  pxt 
psnt(2)  =  pyt 
psnt(3)  =  pzt 
tit  =  volt(l) 
t2t  =  velt(2) 
velt(l)  =  pxt 
velt(2)  =  pyt 
velt(3)  =  pzt 

END  SUB 

FUNCTION  gmstt  (tut,  tmt)   '   03-03-88.     Rev  03-03-88. 

'  Compute  Greenwich  mean  sidereal  time. 

'  tut  =  number  of  centuries  of  36S25  days  of  universal  time  elapsed  since 

'        2000  Jan  1,  12h  UT1  (JD  2451545  UT1). 

'  tmt  =  time  of  day 

'  The  Astronomical  Almanac,  1984,  pp  S13-S15. 

'  Almanac  for  Computers  1988,  pp  B2-B3. 

CONST  pi  =  3.141592653589793t,  rad  =  pi  /  180t 

CONST  cO  =  24110. 54841t,  cl  =  8640184 . 812866t ,  c2  =  9 . 310400000000001D-02 

CONST  c3  =  -.0000062#,  hourpersec  =  It  /  3600t,  meantoapp  =  -.00029t 

CONST  dO  =  125.004452t,  dl  =  -.0529538t,  d2  =  .002071t 

t«  =  (((c3  *  tut  +  c2)  *  tut  +  cl)  *  tut  +  cO)  *  hourpersec  +  tmt'  GMST 

'The  following  two  lines  will  convert  GMST  to  CAST 

'lunarnodet  =  (d2  *  tut  +  dl)  *  tut  +  tOt 

'tt  =  tt  +  meantoapp  *  SIN(lunarnodet  *  rad) 
gmstt  =  dmodt(tt,  24t) 

END  FUNCTION 

SUB  ground. track  (rt,  sat.dct,  sat.ltt,  sat.htt)  '09-14-88.     Rev  09-14-38. 

CONST  ae  =  6378. 14t  '  earth's  equatorial  radius  (km.) 

CONST  fl  =  It  /  298.257t        '  earth's  flattening  factor 

CONST  e.sqr  =  (2t  -  fl)  *  fl    '  earth's  ellipticity  squared 
CONST  f2  =  (It  -  fl)  *  (It  -  fl) 

rkmt  =  rt  *  ae 

pst  =  sat .dct 

el:      cpt  =  COS(pst) 

ret  =  ae  *  SQR((lt  -  e.sqr)  /  (It  -  e.sqr  *  cpt  *  cpt)) 

rcrt  =  ret  /  rkmt 

sat.ltt  =  ATN(TAN(pst)  /  f2)     '  satellite  latitude 

snt  =  SIN(sat.ltt  -  pst) 

hrt  =  SQR(lt  -  rcrt  *  rcrt  *  snt  *  Bnt)  -  rcrt  *  COS(sat.ltt  -  pst) 

pnt  =  sat.dct  -  asin(hrt  *  snt) 

IF  (ABS(pnt  -  pst)  >  .0000001)  THEN  pst  =  pnt:  GOTO  el 
sat.htt  =  hrt  *  rkmt  '  satellite  height 

END  SUB 

'**»***»*«*»»»»*«•*»******«»»*********»*******»«•*»*»»***■»«*»******»*»»*»***»* 

SUB  hr.min.sec  (xt,  hrX,  minX,  sect,  nX) 

'  Convert  decimal  hours  to:     hh  mm  ss      if  nX  =  0 , 

'  hh  mm  ss .f     if  nX  =  1 , 

'  hh  mm  ss.ff    if  nX  =  2, 

'  hh  mm  ss.fff   if  nX  =  3. 

CONST  tO  =  It,  tl  =  lOt,  t2  =  lOOt,  t3  =  lOOOt 
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CONST  rO  =  It  /  7200f,  rl  =  rO  /  tl ,  r2  =  rl  /  t2 ,  r3  =  r2  /  t3 

SELECT  CASE  nX 
CASE  0 

round  =  rO 

trunc  =  tO 
CASE  1 

round  =  rl 

trunc  =  tl 
CASE  2 

round  =  r2 

trunc  =  t2 
CASE  3 

round  =  r3 

trunc  =  t3 
CASE  ELSE 

PRINT 

PRINT  "Error  in  hr. min. sec  subroutine.  Must  have  0  <=  nY.  <=  3." 

PRINT  "nX  =  ";  n% 

PRINT 

EXIT  SUB 
END  SELECT 

vt  =  ABS(xt)  +  round 

hr'/.  =  INT(vf) 

vf  =  (v«  -  hr'/.)  *  60# 

min'/.  =  INT(vf) 

sec*  =  INT((60t  *  (v#  -  min'/.)  *  trunc))  /  trunc 

END  SUB 

'*«********«»*»«»«********»*»***«»*****»«******»»**»»»«.». -•»..•-•••-.•■  ■■••... 
SUB  jd. to. date  (jdf(),  mt ,  dft ,  yft ,  ht ,  mo$,  day$)  '  02-24-88   Rev.  03-02-88 
1  Convert  Julian  Day  Number  to  month  number,  day,  year,  hour,  month  name  and 
'     day  of  the  week.  Limited  to  the  Gregorian  Calendar. 
'  Fliegel  ft  VanFlandern,  Comm.  ACM,  Vol.  11,  No.  10,  Oct.  1968,  pg .  657. 

CONST  w$  =  "MonTueWedThuFriSatSun" 

CONST  dd$  =  "JanFebMarAprMayJunJulAugSepOctNovDec" 

j4ft  =  INT(jdt(l)  +  jdt(2)  ♦  .5) 

ht  =  ((jdt(l)  -  j4ft)  ♦  jdt(2)  +  .5)  *  24t 

It  =  INT(j4ft  +  68569) 

nft  =  INT(4  *  1ft  /  146097) 

lft   =   lft    -    INT((146097   *   nft   +   3)    /  4) 

yft  =  INT(4000  *  (lft  ♦  1)  /  1461001) 

lft  =  It  -  INT(1461  *  yt  /  4)  ♦  31 

mt  =  INT(80  *  It  /  2447) 

dt  =  It  -  INT(2447  *  mt  /  80) 

It  =  INT(mt  /  11) 

mt  =  mt  ♦  2  -  12  *  It 

yt  =  100  *  (nt  -  49)  ♦  yt  ♦  It 

wnX  =  j4t  -  7  *  INT(j4t  /  7)  +  1 
day$  =  MID$(w$,  3  *  wn%  -  2,  3) 
mo$  =  MID$(dd$,  3  *  mt  -  2,  3) 

END  SUB 

'******************************♦***********************************«*»***»»»*** 

SUB  Julian. day .number  (mt,  dt,  yt,  utt,  jdt(),  djt,  tjt) 

'  02-25-88    Rev.  03-23-88 

'  Conversion  of  calendar  date:  mt  =  month,  dt  =  day,  yt  =  4-digit  year, 

'   and  utt  =  Universal  Time,  to  Julian  Day  Number.  Year  must  be  between  1801 

'   and  2099. 
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'  Almanac  for  Computers,  1988,  pg.  B-2. 

Note  that  jd#  is  a  subscripted  double  precision  array.  At  any  time,  the 

Julian  Day  Number  is  given  by  jdf(l)  ♦  jdf(2). 
For  ease  of  programming,  the  user  may  put  the  entire  epoch  in  jdt(l)  and  set 

jdt(2)  =  Of. 
For  maximum  accuracy,  set  jdt(l)  =  the  most  recent  midnight  at  or  before  the 

epoch  and  set  jdf(2)  =  the  fractional  part  of  a  day  elapsed  between  jdf(l) 

and  the  epoch.  As  an  example,  compute  the  number  of  days  from  epoch 

J2000.0  =  JD  2451545.0  using  the  code: 

dayst  =  (jdf(l)  -  2451545. Of)  ♦  jdt(2) 

Note,  particularly,  the  extra  pair  of  parenthesis. 

djf  =  dayB  and  fraction  from  J2000.0.  tjf  =  Julian  centuries  from  J2000.0 

CONST  cl  =  1721013.51,  c2  =  190002. 5f,  jcent  =  36525t 

IF  (1801  <=  yft)  AND  (y*  <=  2099)  THEN 

jd«(l)  =  367  *  yft  -  INT(7  *  (y*  ♦  INT((mft  ♦  9)  /  12))  /  4)  _ 

+  INT((275  *  m*)  /  9)  +  dft  ♦  cl  -  .5t  *  SGN(100  *  yft  +  mft  -  c2)  +  .5* 
jdf(2)  =  utt  /  24 

djf  =  (jdf(l)  -  2451S45t)  ♦  jdt(2) 
tjf  =  dj#  /  jcent 
ELSE 

PRINT 

PRINT  "Error  in  Julian .day .number.  Year  is  not  between  1801  and  2099." 
PRINT  "Year  =  ";  yfc 
PRINT 
END  IF 
END  SUB 

SUB  kepler  (mf,  ecf,  eat)  '  02-25-88.   Rev  02-26-88 

'High  precision,  non-iterative,  solution  to  Kepler's  equation  for  the  ellipse. 
'Accuracy  is  10E-18  using  all  terms  or  10E-15  if  the  last  term  is  omitted. 
'Although  not  as  compact  as  the  Newton-Raphson  procedure,  one  square  root,  one 
'cube  root,  and  only  two  trigonometric  functions  are  evaluated  thus  making  this 
'procedure  extremely  fast.  Can  be  extended  to  the  hyperbolic  case  -  see  Ref . 
'Ref.:  S.  Mikkalo,  Cel.  Mech.,  Vol.  40,  1987,  pp  329-334. 

'mf  =  mean  anomaly,  ecf  =  eccentricity  and  eat  =  eccentric  anomaly. 

CONST  pi  =  3.141592653589793#,  twopi  =  pi  +  pi,  third  =  IS  /  3« 
CONST  r3  =  IS  /  6f ,  r4  =  It  /  24f,  r5  =  It  /  120t 

IF  ect  >=  If  OR  ecf  <  Of  THEN 

PRINT 

PRINT  "Error  in  Kepler.  The  eccentricity  is  not  in  the  [0,1)  interval." 

PRINT  "Eccentricity  =  ";  ecf 

PRINT 
ELSE 

mf  =  mf  -  twopi  *  INT(mf  /  twopi)'  assure  that  mf  is 

IF  mf  >  pi  THEN  mf  =  mf  -  twopi'     in  the  (-pi, pi)  interval 

'Compute  the  initial  approximation: 

gt  =  1  /  (4  *  ecf  ♦  .5) 

af  =  (1  -  ect)  *  gt 

bt  =  .5  *  mf  *  gt 

argt  =  ABS(bt)  ♦  SQR(bf  *  bt  ♦  at  *  af  *  af) 

zf  =  argt  *  third 

IF  bt  <  0  THEN  zt  =  -zf 

st  =  zt  -  af  /  zt 

st  =  st  -  .078  *  st  "  5  /  (1  ♦  ect) 

eat  =  mf  ♦  ecf  *  st  *  (3  -  4  *  st  *  st ) '   abs(dea/ea)  <=  0.002 

36 


'Refine 

f2*  =  e 

f3t  =  e 

fOt  =  e 

lit  =  1 

f4f  =  - 

f5t  =  - 

'ThiB  n 

dt  =  -f 

'Use  mo 

dt  =  -f 

dt  =  -f 

dt  =  -f 

dt  =  -f 

* 

eat  =  e 

END 

IF 

END 

SUB 

the  approximation: 
ct  *  SIN(eat) 
ct  *  COS(eat) 
at  -  f2t  -  mt 
t  -  f3t 
f2t 
f3t 

ext  correction  term  is  sufficient  if  ecc  <  0.25 
Ot  /  fit 

re  of  the  following  terms  to  increase  accuracy 
Ot  /  (fit  ♦  .St   *  f2t  *  dt) 

Ot  /  (fit  ♦  dt  *  (.5i  *  f2t  ♦  dt  *  r3  *  f3t)) 
Ot  /  (fit  ♦  dt  *  (.51  *  f2t  ♦  dt  *  (r3  *  f3t  ♦ 
Ot  /  (fit  ♦  dt  *  (.5t  *  f2t  ♦  dt  *  (r3  *  f3t  ♦ 

(r4  *  f4t  +  dt  *  r5  *  f5t)))) 
at  +  dt 


dt  *  r4  *  f4t))) 
dt  . 


SUB  output4  (kX,  yrft,  mo*,  dak,  time. of. day,  It,  In,  ht ,  el,  az ,  rng,  _ 
langle,  head,  sta.id$,  sat.id$) 

CONST  pi  =  3.141592653589793t,  rad  =  pi  /  180t,  tp  =  pi  ♦  pi 
c$  =  CHR$(34) 


SELECT  CASE  k'/. 
CASE  0 


PRINT  " 

PRINT  "  look" 

PRINT  "year  mo  da  hr  mn   sec    lat 

PRINT  "azim   range  angle   head" 

PRINT 


long 


height 
(km)    elev 


CASE  1 

CALL 
It  = 
In  = 
la  = 
hd  = 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
CASE  ELSE 
PRINT 
PRINT 
PRINT 

END  SELECT 

END  SUB 


r .min . sec(time . of .dayt ,  hrX ,  minX,  sect,  1) 

t  /  rad 

mod(ln  /  rad,  360t) 

angle 

ead  /  rad 

"tttt   tt   tt   tt  tt   tt.t    ";  yrft;  moft ;  dai ; 


USING 

USING 

USING 

*3, 

*3, 

*3, 

«3, 

t3, 

#3, 

#3, 

«3, 


"ttt.tt    ttt.tt    ttttt. t    ttt.tt 
"ttt.tt    ttttt   tt.tt   ttt.tt" 


az 
c$ 
c$ 


It;  In; 

rng;  la 
c$;  sat 

c$;  mot 


hr'/.;  min'/.; 
ht;  el; 
hd 

id$ 
c$, 


c$; 


USING  "!ft!  !A!  ";  c$;  sta.id$ 
USING  "  \tttt\     !**!  ";  c$;  yrft 
USING  "\tt>.     !*S!  ";  c$ ;  daft;  c$;  c$ ;  hrX;  c$ ; 
USING  "\tt>.     \tt.t\    ";  c$;  minX;  c$;  c$;  sec;  c$ ; 
USING  "\ttt.tt\     !•##.##!  ";  c$;  It;  c$;  c$ ;  In;  c$; 
USING  "\ ttttt. t!  \ttt.tt\    ";  c$;  ht ;  c$ ;  c$;  el;  c$; 
USING  "Ittf.tt!  \ttttt\    ";  c$;  az ;  c$;  c$;  rng;  c$; 
USING  ••'.tt.tt*     \ttt.tt\";    c$;  la;  c$;  c$;  hd;  c$ 


"Error  in  output4  subroutine.  Must  have  0  <=  k%  <=  1." 
"ky.  =  ";  kX 


SUB  pos. update  (iy.,  elemt(),  timet,  sat.xt(),  sat.xdt(),  elemut(),  rt ,  eat)  _ 

STATIC 
'  09-07-88.  Rev.  03-09-89. 

'  Case  1 :  Initialize  orbital  elements  for  time  of  epoch. 

'  Case  2:  Update  orbital  elements  and  compute  in-plane  position  ft  velocity. 

'  Simplest  of  NAVSPASUR  methods.  Accurate  to  10-15  n.mi.  within  20  days 
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of  epoch 

'     elem(l) 

= 

amajor 

'      elem(2) 

= 

ecc 

'      elem(3) 

= 

tzero 

'      elem(4) 

= 

aincl 

'      elem(5) 

= 

anode 

'      elem(6) 

= 

peri 

'      elem(7) 

= 

motion 

'      elem(8) 

= 

manomaly 
Hnrav 

semi -major  axis 

eccentricity 

time  of  perifocal  passage 

inclination  (radians) 

longitude  of  ascending  node  (radians) 

argument  of  perifocus  (radians) 

mean  motion  (revolutions/day) 

mean  anomaly  (radians) 

1st  decay  constant  (radians/minute~2) 

CONST  pi  =  3.1415926535897931,  rad  =  pi  /  180#,  tsopi  *  pi  +  pi 

CONST  c  =  .00081197385t   '   0.75*j2 

CONST  tuothird  =  It   I   39 

CONST  ae  =  6378.141  '  earth's  equatorial  radius  (km.) 

CONST  GE  =  398600. 5#  '  geo .  grav.  const.  (km)-3/(sec)*2 

X  =  60f  *  SQR(GE  /  ae)  /  ae' (eru)"(3/2)/min 
CONST  we  =  1.0027379093507951  *  twopi  /  14401 'sidereal  period  in  radians/min 

SELECT  CASE  1% 
CASE  1 

elem#(l)  =  (k  /  elemf(7))  *  twothird 

a2f  =  elemf(l)  "  2 

ecc2#  =  If  -  elem«(2)  *  2 

ci#  =  C0S(elemf (4)) 

c2if  =  cif  *  cit 

elemt(l)  =  elemt(l)  *  (If  +  c  *  (3t  *  c2it  -  It)  /  _ 
(a2f  *  ecc2t  "  1.5t))  *  twothird 

elem#(3)  =  timet  -  elemf(8)  /  elemt(7) 

aecc22t  =  (elemt(l)  *  ecc2t)  *  2 

peridotf  =  c  *  (5t  *  c2if  -  It)  /  aecc22t 

nodedotf  =  -2f  *  c  *  cit  /  aecc22t 

adotf  =  -2t  *  twothird  *  elemf(l)  *  m2t  /  elemf(7) 

FOR  i%  =  1  TO  9 

elemuf(iy.)  =  elemf(i*/,) 

NEXT 
CASE  2 

ma.mOf  =  (elemf(9)  *  timet  +  elemf(7))  *  timet 

mat  =  ma.mOf  ♦  elemf(8) 

elemut(6)  =  elemt(6)  +  peridotf  *  ma.mOf 

elemut(5)  =  elemt(5)  +  nodedotf  *  ma.mOf  -  we  *  timet 

elemuf(l)  =  elemf(l)  ♦  adotf  *  timet 

CALL  kepler(maf,  elemf(2),  eat) 
sOf  =  SIN(eaf) 
cOf  =  COS(eaf) 

elemuf(8)  =  mat 

elemuf(3)  =  timet  -  mat  /  elemf(7) 

'  in-plane  position  and  velocity 

sat.xt(l)  =  elemuf(l)  *  (cOt  -  elemt(2)) 

sat.xf(2)  =  elemut(l)  *  SQR(ecc2t)  *  sOt 

sat.xt(3)  =  Of 

rf  =  elemut(l)  *  (It  -  elemt(2)  *  cOt) 

edt  =  elemut(l)  *  elemt(7)  /  rt 

sat.xdt(l)  =  -elemut(l)  *  edt  *  sOf 

sat.xdf(2)  =  elemuf(l)  *  edt  *  SQR(ecc2t)  *  cOt 

sat.xdt(3)  =  Of 

'  equatorial  rectangular  coordinates 

CALL  euler(elemu(6) ,  elemu(5),  elem(4) ,  sat.x(),  sat.xdO) 

'  correct  velocity  for  earth's  rotation 
sat.xdt(l)  =  sat.xdt(l)  ♦  we  *  sat.x(2) 
sat.xdf(2)  =  8at.xdf(2)  -  we  *  sat.x(l) 
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CASE  ELSE 

PRINT 

PRINT  "Error  in  pos  . update.  First  argument  must  be  1  or  2  only." 

PRINT  "Argument  =  ";  1% 

PRINT 
END  SELECT 
END  SUB 

'*******************»**********♦»*»**»********»#*»**»»»»*»»*»»»»»#*»»»»»**»*«*. 
SUB  rise. set  (zv#(),  sta.RtQ,  set,  tOS,  eat,  elemutO) 

Compute  eccentric  anomaly  for  the  time  of  riBe  or  set  of  a  satellite 
02-04-89.  Rev.  02-08-89  C  1700. 

NOTE:  To  use  this  routine,  the  time  of  culmination  MUST  be  computed  first, 
The  entering  eat  argument  should  be  ea.cult  -/♦  0.2  radians  (about  11 
degrees)  for  rise/set. 

See  the  introductory  comments  in  SUB  culminate  for  nomenclature  and 
notation . 

CONST  pi  =  3.1415926S3589793t,  rad  =  pi  /  180t,  twopi  b  pi  +  pi 

amajt  =  elemuf(l) 
ecct  =  elemut(2) 
captf  =  elemu#(3) 
inclt  =  elemut(4) 
nodet  =  elemut(5) 
peril  =  elemu#(6) 
motf  =  elemuf(7) 

DIM  pvt(3),  pvdt(3),  qvt(3).  qvdf(3),  rvecdf(3),  rhof(3) 

coswf  =  C0S(peri#) 
sinwt  =  SIN(perit) 
cosit  =  COS(inclf) 
sinit  =  SIN(mclt) 


corrf 
DO 


coseaf  =  COS(eat) 
ecoseaf  =  eccf  *  coseat 
sineat  =  SIN(eaf) 
esineat  =  eccl  *  sineat 

gammaf  =  nodet  -  wet  *  ((eat  -  esineat)  /  mott  +  captt  -  tOt) 
gammadt  =  -wet  *  (It  -  ecoseat)  /  mott 

cosgt  =  COS(gammat) 
smgt  s  SIN(gammat) 

rtecct  =  amajt  *  SQR(lt  -  ecct  *  ecct) 

xwt  a  amajt  *  (coseat  -  ecct) 
xwdt  =  -amajt  *  sineat 
ywt  =  rtecct  *  sineat 
ywdt  =  rtecct  *  coseat 

pvt(l)  =  coswt  *  cosgt  -  sinwt  *  singt  *  cosit 
pvt(2)  -  coswt  *  singt  ♦  sinwt  *  cosgt  *  cosit 
pvt(3)  =  sinwt  *  sinit 

qvt(l)  =  -sinwt  *  cosgt  -  coswt  *  singt  *  cosit 
qvt(2)  =  -sinwt  *  singt  ♦  coswt  *  cosgt  *  cosit 
qvt(3)  ~    coswt  *  sinit 


39 


pvdf(l)  =  -pvf(2)  *  gammadf 

pvdf(2)  =  pvt(l)  *  gammadf 

pvdt(3)  =  Of 

qvdf(l)  =  -qv#(2)  *  gammadf 

qvdf(2)  =  qvf(l)  *  gammadf 

qvdf(3)  =  Of 

FOR  i%   -    1  TO  3 

rhof(iX)  =  xwf  *  pvf(iX)  +  ywf  *  qvf(iX)  -  sta.Rf(iy.) 
rvecdf(iX)  =  xwdf  *  pvf(iX)  ♦  xwf  *  pvdf(iX)  ♦  ywdf  *  qvf(i%)  _ 
♦  ywf  *  qvdf(iX) 

NEXT 

fef  =  dotf(zvf(),  rhofO) 
fedf  =  dotf(zvf(),  rvecdfO) 

corrf  =  fef  /  fedf 
eaf  =  eaf  -  corrf 

LOOP  UNTIL  ABS(corrf)  <  .OOOOlf 
END  SUB 

>»*****.*»***»»»***»*************************»****»«*»*«***»*♦***»»»»*«*****»** 
SUB  sat. head  (rf,  sat.xf(),  sat.xdf(),  headf)  '   01-18-89.  Rev.  03-12-89. 
'satellite  heading 

CONST  pi  =  3.141592653589793f ,  twopi  =  pi  +  pi 

sheadf  =  (sat.xf(l)  *  sat.xdf(2)  -  sat.xf(2)  *  sat.xdf(l))  /  rf 

cheadf  =  sat.xdf(3)  -  (sat.xf(3)  *  dotf (sat .xf () ,  sat.xdf()))  /  rf  "  2 
headf  =  atan2f (sheadf ,  cheadf) 

IF  headf  <  Of  THEN  headf  =  headf  +  twopi 

END  SUB 

SUB  sat . sphere. coord  (sat.x(),  r,  t j ,  tm,  sat.dc,  sat. In,  sat.ra) 
'compute  spherical  coordinates 

CONST  pi  =  3.141592653589793f ,  rad  =  pi  /  180f,  twopi  =  pi  ♦  pi 

sat.dc  =  asin(sat .x(3)  /  r)  '  satellite  declination 

sat. In  =  atan2(sat .x(2) ,  sat.x(l))  '  satellite  longitude 

ru  =  15f  *  gmst(tj,  tm)  *  rad  '  Greenwich  mean  sidereal  time 

sat.ra  =  dmod(sat.ln  +  ru ,  twopi)  '  satellite  right  ascension 

END  SUB 

SUB  sta.p.coord  (ltf,  lnf ,  htf,  xf())  '   08-06-88.    Rev  02-04-89. 
'  convert  station,  lat,  long  t   height  (ft.)  from  geographic 
'    to  geocentric  rectangular  -  x(.) 

CONST  nmi.ft  =  .3048f  /  1852f  '  n. mi. /foot 

CONST  km. ft  =  .3048f  /  lOOOf  '  km. /foot 

CONST  ae  =  6378. 14f  '  earth's  equatorial  radius  (km.) 

CONST  fl  =  If  /  298.257f  '  earth's  flattening  factor 

CONST  e.sqr  =  (2f  -  fl)  *  fl  '  earth's  ellipticity  squared 

ht.eruf  =  htf  *  km. ft  /  ae 

sltf  =  SIN(ltf) 

cltf  =  COS(ltf) 

slnf  =  SIN(lnf) 

clnf  =  COS(lnf) 

nf  =  If  /  SQR(lf  -  e.sqr  *  sit  *  sit) 

glf  =  nf  ♦  ht.eruf 
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g2t  ■  nf  *  (It  -  e.sqr)  ♦  ht . eru« 

'radius  vector  for  flattened  earth 
x#(l,  1)  =  glf  *  cltt  *  clnt 
xt(2,  1)  =  gl«  *  cltt  *  slnt 
xt(3,  1)  =  g2t  *  sltt 

'Spherical  earth  unit  vectors: 

'1)  local  perpendicular  vector 
xt(l,  2)  =  cltt  *  clnt 
xt(2,  2)  =  cltt  *  slnt 
xt(3,  2)  =  sltt 

'2)  east  vector 
xt(l,  3)  =  -slnt 
xt(2,  3)  =  clnt 
xt(3,  3)  =  Ot 

'3)  north  vector 
xt(l,  4)  =  -sltt  *  clnt 
xt(2,  4)  =  -sltt  *  slnt 
xt(3,  4)  =  cltt 

END  SUB 

SUB  sta.sat  (sta.xt(),  sat.xt(),  rngt ,  azt,  elt,  langlet) 

08-06-88.  Rev.  08-31-88. 

'compute  station-satellite  relations.  Range,  azimuth,  elevation  and  satellite 
'  "look-angle"  (angle  between  satellite  sub  point  and  station) . 

CONST  pi  =  3.141592653589793t,  rad  =  pi  /  180t,  twopi  =  pi  ♦  pi 

*rho()  =  sat  coords  minus  sta  coords,  rhodots  are  the  rho  components  projected 
'   in  the  station  local  north-east-radial  system, 
DIM  rhot(3),  rhodotst(3) 

rngt  =  Ot 

FOR  i'/.  =  1  TO  3 

rt  =  sat.xt(i*/.)  -  sta.xtd?.,  1) 

rhot(i'/.)  =  rt 

rngt  =  rngt  +  rt  *  rt 
NEXT 

rngt  =  SQR(rngt) 
FOR  i%   =  1  TO  3 

rhot(iy.)  =  rhot(i'/.)  /  rngt 
NEXT 

FOR  iX  =  1  TO  3 

at  =  Ot 

FOR  jX  =  1  TO  3 

at  =  at  +  rhot(jy.)  *  sta.xt(jy,,  if.   ♦  1) 

NEXT 

rhodotst(iy.)  =  at 
NEXT 

azt  =  atan2t(rhodotst(2) ,  rhodotst(3))  /  rad 
IF  azt  <  Ot  THEN  azt  =  azt  «•  360t 

at  =  rhodotst(l) 

elt  =  ATN(at  /  SQR(lt  -  at  *  at))  /  rad 

at  =  Ot 
bt  =  Ot 
FOR  it   =    i  TO  3 
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ct  =  sat.x#(iy.) 

at  =  a*  +  rhot(iX)  *  ct 

bt  =  bt  ♦  ct  *  ct 

NEXT 

langlet  =  acost(at  /  SqR(bt))  /  rad 

END  SUB 

SUB  time. update  (k%,  delta,  tt,  tmt,  jdt,  jdt(),  tjt,  time,  of  .day ,  mJt ,  dt ,  _ 
yft,  ht,  mo$,  day$) 

CONST  twentyfour  =  24t,  sixty  =  60t,  Jul. cent  =  36525t 

tt  =  delta. tt  'delta. tt  in  hours  if  kX  =  0,  tt  in  hours. 

IF  kX  >  0  THEN  tt  =  tt  /  sixty    'delta. tt  in  minutes  if  kX  >  0,  tt  in  hours. 

tmt  =  tmt  ♦  tt 

jdt(2)  =  jdt(2)  ♦  tt  /  twentyfour 

jdt  =  jdt(l)  ♦  jdt(2) 

tjt  =  tjt  ♦  tt  /  twentyfour  /  jul.cent 

time. of. day  =  time. of. day  +  tt 

IF  time. of .day  >=  twentyfour  THEN 

time. of. day  =  time. of. day  -  twentyfour 

jdt(l)  =  jdt(l)  +  It:  jdt(2)  =  jdt(2)  -  It 

CALL  jd.to.date(jdt() ,  mfc,  dft,  yft ,  ht,  mo$,  day$) 
END  IF 
IF  time. of. day  <  Ot  THEN 

time. of .day  =  time. of .day  +  twentyfour 

jd#(l)  =  jdt(l)  -  It:  jdt(2)  =  jdt(2)  +  It 

CALL  jd.to.date(jdt() ,  mft,  dft,  yt,  ht,  mo$ ,  day$) 
END  IF 

END  SUB 

SUB  yrmoday  (v$,  yrft,  moft ,  dyft)  '   03-21-88.     Rev  03-21-88. 
'  Convert  yyyy.mmdd  to  yr&=yyyy ,  moft=mm,  and  dyft=dd 

IF  INSTR(v$,  ".")  <>  5  OR  LEN(v$)  <>  9  THEN 

PRINT 

PRINT  "Error  in  yrmoday.  v$  must  be  of  the  form  yyyy.mmdd" 

PRINT  "v$  =  ";  v$ 

PRINT 
ELSE 

yr*  =  VAL(LEFT$(v$,  4)) 

mo*  =  VAL(MID$(v$,  6,  2)) 

dyft  =  VAL(RIGHT$(v$,  2)) 
END  IF 
END  SUB 
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APPENDIX  E.  THE  HIGH  ACCURACY  KINEMATICS  ROUTINE. 

This  appendix  contains  a  listing  of  the  alternate  pos. update  routine.  It  is  a  BASIC  ver- 
sion of  the  highly  accurate  PPT2  subroutine  provided  by  NAVSPASUR  [Ref.  1]  with  their 
SHOWALL  program.  It  is  believed  that  it  has  been  accurately  translated,  however  any  errors 
are  solely  those  of  the  author.  Known  discrepacies  are  the  difference  in  the  WGS-72  con- 
stants used  in  SHOWALL  and  the  IAU  1976  constants  used  here  [see  section  III.  SOURCES, 
in  this  report  for  details]. 

»***»*»»»««***«»»»»»»*»»»•»»*«»»»*•*»»»»•»««»»«»••»«••••»••• ............ 

SUB  pos. update  (ind%,  elemf(),  timet,  sat.x*(),  sat.xdi(),  elemuf(), 
r*.  eat)  STATIC    Rev.  09-26-89  0  0900 

Case  1:  Initialize  orbital  elements  for  time  of  epoch. 

Case  2:  Update  orbital  elements  and  compute  in-plane  position  b   velocity. 

Comments  from  original  FORTRAN  version: 
implicit  real*8(a-h ,o-z) 

model  reference-Astronomical  Journal  64, No  1274 ,Nov . , ' 59 
solution  of  the  problem  of  artificial  satellite  theory 
with  out  drag... Dirk  Brouwer. 

840329   confirmed  use  of  w=u  x  v   vice   w=u  x  vel 
840329   inserted  cdh.sdh,  cdt.sdt  re  osc  vice  mean  els  at  t 
for  use  only  with  dcext2 ,dcl ,dcia 

frame  for  u.v.w  is  the  rotated  inertial  w/o  coriolis 
830127   new  partials  re  (a ,x ,y ,z ,p ,q) 
830127   improved  (polynomial)  model  for  singularity  in  Brouwer 


*** 
*** 


DIM  a(25) 

DIM  f(25),  osc(10),  kf(10),  cf(10),  bs(3,  4),  u(3)  ,  v(3)  ,  w(3),  vel(3) 

DIM  ar(6,  8),  br(3,  6),  cr(3,  6),  dv(3) 


'  *** 
'  *** 
'  *** 
'  *** 
'  *** 


The  commented  beta  value  is  the  WGS-72  value .however  using  this 

value,  one  obtains  a  value  for  a(9)  different  then  what  we  had 

been  using,  thus  causing  a  problem  with  time . therefore  to  retain 

previous  value  of  a(9),one  solves  for  beta  using  the  old  value 

of  a(9)  . 

betal=398600.5d0 

betal  =  398597.625795881,  erkm  =  6378. 135#,  flat  =  298. 26t 

k  =  .0743669161# 


'  IAU  1976: 

CONST  twopi  =  6.283185307179586t 
CONST  we  =  1.0027379093507951  *  twopi 
CONST  ae  =  6378. 14* ,  erkm  =  ae 
CONST  betal  =  398600.5* 

k  =  60*  *  SQR(betal  /  ae)  /  ae 
CONST  flat  =  1*  /  298.257* 


1440*' sidereal  period  in  radians/mm 

'  earth's  equatorial  radius  (km.) 

'  geo.  grav.  const,  (km) *3/(sec) "2 

'(eru)*(3/2)/min 

'  earth's  flattening  factor 


'     k-terms      [Note:  These  may  or  may  not  agree  with  IAU  1976] 
CONST  c20  =  -.0004841605*,  c30  =  .00000095958* 
CONST  c40  =  .00000055199*.  c50  =  .000000065875* 


CONST  twotrd  =  .666666666667*.  fortrd  = 
CONST  tentrd  =  3.333333333333*.  ar  =  48 


1.333333333333* 
*  0* 


'IF  indY.  <>  1  THEN  GOTO  secondcase 
SELECT  CASE  ind% 
CASE  1 

kZy.  =  i 

a(l)  =  0*'   [Note:  =  85-01-01  in  SHOWALL.  Set  to  zero  here.] 
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a(3)  =  -.59   *   c20  *  SQR(5t) 

a(4)  =  c30  *  SqR(7t) 

a(5)  =  .375  *  c40  *  SQR(9t) 

a(6)  =  c50  *  SQR(llt) 

'      tso  pi 

a(7)  =  tHopi 

a(l6)  =  If  /  a(7) 

'      hergs/day  ,  secs/herg,  mins/herg 

a(9)  =  erkm  *  SqR(erkm  /  betal) 

a(8)  =  86400t  /  a(9) 

a(l7)  =  1440f  /  a(8) 

'     we(rad/day),  we-  2  pi,  we(rad/herg) 

a(ll)  =  6.3003880987t 

a(10)  =  a(ll)  -  a(7) 

a(l2)  =  a(ll)  /  a(8) 

'      earth  flattening 

a(13)  =  (2t  -  It  /  flat)  /  flat 

'     sm/er,  km/er,  nm/er 

a(20)  =  erkm  /  1. 6093441 

a(21)  =  erkm 

a(22)  =  erkm  /  1.852# 

'      deg/rad 

a(23)  =  360t  /  twopi 

'      range  rate/er/herg  to  cycles/second  -  conversion 

a(24)  =  a(21)  *  216980000#  /  (a(9)  *  299792. 5#) 

'      fence  plane  displac  from  earth  center 

a(25)  =  .00311 

'  Patch  to  interface  PPT2  and  SATSTA  variables: 

f(l)  =  elem(8)  'mean  anomaly 

f(2)  =  elem(7)  *  a(17)  'mean  motion 

f(3)  =  elem(9)  *  a(l7)  *  2'first  decay  constant 

f(4)  =  Of  'second  decay  constant 

f(5)  =  elem(2)  'eccentricity 

f(6)  =  elem(6)  'arg  of  perigee 

f(7)  =  elem(5)  'node 

f(8)  =  C0S(elem(4))     'cosine  inclination 

f(9)  =  Of     '  epoch  in  PPT2,  time  from  epoch  here  (?  but  check  ?) 

'   Used  only  once  to  "recover"  the  remaining  elements 

hi  =  Of 

'***  f(8)  is  cosine  of  inclination 

t9  =  f(8)  *  f(8) 

tl  =  If  -  5t  *  t9 

t2=(1.0d0-exp(-100.0d0*tl"2))/tl 
'     revised  handling  of  divisor  per  navspasur  o2t  memo  nov  '83 
beta  =  lOOt  /  2t  *  11 
p3  =  beta  *  tl  "  2 
pi  =  EXP(-p3) 
p2  =  It  ♦  pi 
plex  =  pi  "  2 
p4  =  it 

p3ex  =  -.5t  *  p3 
FOR  nX  =  2  TO  13 

IF  n%  <=  11  THEN  p2  =  p2  *  (It  +  plex) 

plex  =  plex  *   plex 

p4  =  p4  ♦  p3ex 

p3ex  =  -p3  *  p3ex  /  (nX  +  1) 
NEXT 

t2  =  p2  *  p4  *  beta  *  tl 
cio2  =  SQR(.5t  ♦  .5f  *  f(8)) 
sio2  =  SQR(.5t  -  .5t  *  f(8)) 
tl2  =  3t  *  t9  -  It 
tl3  =  t9  *  t9 
tl4  =  t2  *  tl3 
tl5  =  It  -  t9 
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>***   f(5)  is  eccentricity 

eBqO  =  f(5)  *  f(5) 

eta20  =  1*  -  esqO 

etaO  =  SQR(eta20) 

t6  =  etaO  *  eta20 

t7  =  eta20  *  eta20 

t8  =  etaO  +  If  /  (It  ♦  etaO) 

f(13)  =  f(2)  *  (-twotrd) 

elem(l)  =  f(l3) 

'***  f(2)  is  basically  mean  motion  but  also  contains  the 

'***   secular  term  in  mean  anomaly,  therefore  f(13)  is  (almost)  the 

'***   semi  major  axis. 

>***   This  loop  is  to  remove  the  contribution  of  the  additional 

'***   term  in  mean  motion  for  improving  values  of  semi -major  axis 

>***   and  seculars. 

>***   a(3)  to  a(6)  are  'k'  terms  wgs  72 

FOR  iX  =  1  TO  5 

fsl3  =  f(13)  *  f(13)  *  t7 

fsll  =  a(3)  /  fsl3 

fsl3  =  a(5)  /  (fsl3  *  fsl3) 

tt  =  1*  +  1.5t  *  fsll  *  etaO  *  tl2 

uu  =  -I5t  ♦  16f  *  etaO  ♦  25t  *  eta20 

uu  =  uu  +  (308  -  96f  *  etaO  -  90t  *  eta20)  *  t9 

uu  =  uu  +  (105*  +  144t  *  etaO  ♦  25t  *  eta20)  *  tl3 

tt  =  tt  +  .09375f  *  fsll  "  2  *  etaO  *  uu 

tt  *  tt  ♦  .9375f  *  fsl3  *  etaO  *  esqO  *  (3t  -  30f  *  t9  ♦  35t  *  tl3) 

f(l3)  =  (tt  /  f(2))  "  twotrd 
NEXT 

'***   f(ll)  will  be  rate  of  change  of  argument  of  perigee 
uu  =  -35f  +  24f  *  etaO  +  25t  *  eta20 
uu  =  uu  +  (90t  -  192f  *  etaO  -  126t  *  eta20)  *  t9 
uu  =  uu  ♦  (385t  +  360t  *  etaO  +  45f  *  eta20)  *  tl3 
tt  s  -1.5t  *  fsll  *  tl  +  .09375f  *  fsll  *  2  *  uu 
uu  =  21*  -  9«  *  eta20  +  (-270t  +  126#  *  eta20)  *  t9 
uu  =  uu  ♦  (385*  -  189*  *  eta20)  *  tl3 
f(ll)  =  tt  +  .3125t  *  fsl3  *  uu 

'***   f(12)  will  be  rate  of  change  of  right  ascention 
uu  =  -5*  +  12«  *  etaO  ♦  9*  *  eta20 
uu  =  uu  -  (35*  +  36#  *  etaO  +  5*  *  eta20)  *  t9 
tt  =  -3#  *  fsll  *  f(8)  ♦  .375«  *  fsll  ■  2  *  f(8)  *  uu 
f(12)  =  tt  ♦  1.25*  *  fsl3  *  f(8)  *  (5*  -  3*  *  eta20)  *  (3*  -  7*  -  t9) 
'***   f(15)  =  sme  inclination.f (14)is  a  dot,f(16)  is  eccen.dot 
f(15)  =  SQR(tl5) 

f(l4)  =  -fortrd  *  f(3)  *  1(13)  /  f(2) 
f(16)  =  f(14)  *  f(5)  *  eta20  /  f(13) 
ql  =  f(2)  *  f(13)  "  1.5* 
f(U)  =  f(ll)  /  ql 
f(12)  =  f(12)  /  al 

fsl2  =  a(4)  /  (a(3)  *  f(13)  *  eta20) 
fsl3  =  fsl3  /  fsll  *  tentrd 

fsl4  =  a(6)  /  (a(3)  *  f(13)  "  3  *  eta20  *  t7) 
ql  =  .125*  *  (fsll  *  (1*  -  11*  *  t9  -  40*  *  tl4)  -  fsl3  _ 

*  (1*  -  3*  *  t9  -  8*  *  tl4)) 
p5  =  1*  ♦  t2  *  (8*  *  t9  ♦  20*  *    tl4) 
p2  =  1*  +  2*  *  p5 

q2  =  .125*  *  esqO  *  f(8)  *  (fsll  *  (1*  ♦  10*  *  p5)  -  fsl3  *  p2) 
p2  =  .46875*  *  p2  *  f(5)  *   f(8)  *  f(15)  *  (4*  ♦  3*  *  esqO)  *  fsl4 
p3  =  fsl4  *  (1*  -  9*  *  t9  -  24*  »  tl4) 
q5  =  .25*  *  (fsl2  ♦  .3125*  *  (4*  ♦  3*  *  esqO)  *  p3) 
p3  =  .15625*  *  f(5)  *  f(l5)  *  p3 

p4  =  .030381944*  *  f(5)  *  fsl4  *  (1*  -  5*  *  t9  -  16*  *  tl4) 
p5  =  .060763889*  *  f(5)  *  esqO  «  f(8)  *  f(15)  *  fsl4  *  (1*  +  4*  *  p5) 
vlel  =  f(5)  *  eta20  *  ql 
vlhli  =  -f(15)  *  q2 
tt  =  (t6  -  1*)  *  ql  -  q2  ♦  25*  *  esqO  *  tl4  *  t9  *  t2  *  (fsll  -  .2*  *  fsl3) 
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vlsl 

vle2 

vlh2i 

vls2 

vle3 

vlh3i 

vls3 

vll2 
i 

sinhO 
coshO 
sintO 
costO 

> 

dadm  = 

secul 
dgddm 
dgdde 
dgddi 
dhddra 
dhdde 
dhddi 
dtddx 
dtddy 
dtddp 
dtddq 
dhddx 
dhddy 
dhddp 
dhddq 
> 

daddml 
daddm2 
dedde 
deddx 
deddy 
deddml 
deddm2 
FOR  iV. 
e 
NEXT 


=  tt  -  .0625f  *  esqO  *  (fall  *  (It  -  33t  *  t9  -  200t  *  tl4)  _ 

-  Isl3  *  (It  -  9t  *  t9  -  40t  *  tl4)) 
=  eta20  *  1(15)  *  q5 
=  1(5)  *  1(8)  *  q5  ♦  f(15)  *  p2 
=  f(5)  *  1(15)  *  (t8  ♦  1(8)  /  (It  ♦  f(8)))  *  q5  . 

♦  (lit  ♦  3t  *  esqO  -  3t  *  t6)  *  p3  ♦  (It  -  f(8))  *  p2 
=  -3t  *  1(5)  *  eta20  *  1(15)  *  p4 
=  -eBqO  *  1(8)  *  p4  -  1(15)  *  p5 

t   1(15)  *  (3t  *  t6  -  3t  -  2t  •  esqO  -  esqO  *  1(8)  /  (It  ♦  1(8)))  *  p4 
-  (It  -  1(8))  *  P5 

■  vle2  ♦  3t  *  1(5)  *  eta20  *  p3 
precomps  lor  the  partials 

1(6)  is  argument  of  perigee, 1(7)  is  right  ascension 
=  SIN(1(7)) 
=  C0S(1(7)) 
=  SIN(1(6)  ♦  1(7)) 
=  C0S(1(6)  ♦  1(7)) 
kepler's  law 

twotrd  *  1(13)  /  1(2) 
variations  thru  the  seculars 
=  a(3)  /  (1(13)  *  eta20)  *  2 
=  3t  *  secul  *  tl  /  1(13)  *  dadm 
=  -6t  *  secul  *  tl  *  1(5)  /  eta20 
=  -15t  *  secul  *  1(8)  *  1(15) 
=  6t  *  secul  *  1(8)  /  1(13)  *  dadm 
=  -12t  *  secul  *  1(8)  *  1(5)  /  eta20 
=  3t  *  secul  *  1(15) 

■  (dgdde  +  dhdde)  *  costO 


sintO 

2t  *  coshO  /  cio2 

2t  *  sinhO  /  cio2 


=  (dgdde  ♦  dhdde) 

=  (dgddi  +  dhddi) 

=  (dgddi  +  dhddi) 

■  dhdde  *  costO 

=  dhdde  *  sintO 

=  dhddi  *  2t  *  coshO  /  cio2 

=  dhddi  *  2t  *  sinhO  /  cio2 

variations  thru  m2 

=  lortrd  *  1(3)  *  (1(13)  /  1(2)  -  dadm)  /  1(2) 

=  -fortrd  *  1(13)  /  1(2) 
=  -fortrd  *  (It  -  3t  *  esqO)  *  1(3)  /  1(2) 
=  dedde  *  costO 
=  dedde  *  sintO 

=  lortrd  *  f(5)  *  eta20  *  1(3)  /  1(2)  "  2 

=  -fortrd  *  1(5)  *  eta20  /  1(2) 
1  TO  9 
lemutd'/.)  =  elemt(iy.) 


'  secondcase : 
CASE  ELSE 

tm  =  timet  /  a(l7)    'Convert  minutes  to  Hergs 

hi  =  tm  -  1(9) 

>***  compute  position  as  a  lur.ction  ol  time  (in  call  line) 

osc(9)  =  hi  *  (1(2)  ♦  hi  *  (1(3)  ♦  hi  *  1(4))) 

osc(3)  =  dmod(Kl)  +  osc(9),  twopi) 

tt  =  1(5)  +  1(16)  *  hi 

IF  tt  <  Ot  THEN  tt  =  Ot 

IF  tt  >  .99999t  THEN  tt  =  .99999t 

osc(4)  =  tt 

esq  =  osc(4)  *  2 

eta2  =  It  -  esq 

eta  =  SQR(eta2) 

CALL  kepler(o8c(3) ,  osc(4),  c3) 

osc(l)  =  C0S(c3) 

cl  =  It  -  osc(4)  *  osc(l) 

osc(l)  =  (osc(l)  -  osc(4))  /  cl 
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osc(2)  =  eta  *  SIN(c3)  /  cl 

t84  =  It  ♦  obc(4)  *  obc(1) 

cf3  =  Of 

IF  kzX  <>  0  THEN  cf3  =  a(12) 

tt  =  f(13)  ♦  f(14)  *  hi 

IF  tt  <  It  THEN  tt  =  It 

osc(8)  =  tt 

osc(7)  =  f(8) 

f8i  =  f(l5) 

osc(6)  =  dmod(f(7)  ♦  f(12)  *  osc(9),  tnopi) 

osc(5)  =  dmod(f(6)  ♦  f(ll)  *  osc(9),  tnopi) 

r  =  08c(8)  *  eta2  /  ts4 

osc(6)  =  obc(6)  -  cf3  *  (tm  -  a(l)) 

wl  =  C0S(obc(5)) 

w2  =  SIN(osc(5)) 

w7  =  C0S(osc(6)) 

w8  =  SIN(o8c(6)) 

agda  =  Ot 

agde  =  Ot 

agdi  =  Ot 

agdl  =  Ot 

agdg  =  Ot 

agdh  =  Ot 

w3  =  Hi  *  wl  -  w2  *  w2 

w4  =  2t  *  wl  *  w2 

w5  =  h3  *  wl  -  w4  *  w2 

w6  =  w4  *  wl  +  w3  *  w2 

h9  =  osc(l)  *  osc(l)  -  osc(2)  *  osc(2) 

HlO  =  2t  *  osc(l)  *  osc(2) 

wll  =  w3  *  osc(l)  -  w4  *  osc(2) 

wl2  =  h4  *  osc(l)  +  w3  *  osc(2) 

wl3  =  w3  *  h9  -  h4  *  HlO 

h14  =  h4  *  h9  +  h3  *  HlO 

wl5  =  h13  *  osc(l)  -  wl4  *  08c(2) 

wl6  =  h14  *  osc(l)  +  Hl3  *  osc(2) 

h17  =  atan2(osc(2) ,  osc(l))  ♦  osc(4)  *  osc(2)  -  dmod(osc(3),  twopi) 

wl8  =  C0S(osc(3)) 

wl9  =  SIN(osc(3)) 

w20  =  osc(l)  *  (3t  ♦  osc(4)  *  osc(l)  *  (3t  ♦  osc(4)  *  osc(l))) 

h21  =  3t  *  h14  +  3t  *  osc(4)  *  h12  +  osc(4)  *  h16 

h22  =  ts4  *  (It  ♦  ts4)  /  eta20 

osc(8)  =  osc(8)  *  (It  ♦  fell  /  eta20  *  (tl2  *  (ts4  "  3  -  t6)  _ 

+  3t  *  tl5  *  ts4  "  3  *  h13))  +  agda 
de  =  vlel  *  h3  ♦  vle2  *  h2  ♦  vlo3  *  w6 
di  =  sio2  +  .5t  *  cio2  *  (.5t  *  fall  *  f(8)  *  f(15)  *  (3t  *  wl3  +  3t  *  f(5)  _ 

*  wll  ♦  f(5)  *  h15)  -  de  *  f(5)  *  f(8)  /  f(l5)  /  eta20)  ♦  .St   *   cio2  *  agdi 
de  =  osc(4)  ♦  .St  *  fall  *  (tl2  *  (h20  ♦  f(5)  *  t8)  ♦  3t  *  tl5  *  (h20  ♦  f(5))  _ 

*  h13  -  eta20  *  tl5  *  (3t  *  wll  ♦  h15))  ♦  de  ♦  agde 

dh  =  .5t  /  cio2  *  (-.5t  *  fall  *  f(8)  *  f(15)  *  (6t  *  h17  -  w2l)  ♦  vlhli  *  h4  . 

♦  vlh2i  *  Hi  ♦  vlh3i  *  h5)  ♦  .5t  *  agdh  /  cio2 
dl  =  -.25t  *  t6  *  fell  *  (2t  *  tl2  *  (h22  ♦  It)  *  osc(2)  ♦  3t  *  tl5  *  ((it  _ 

-  h22)  *  h12  ♦  (.333333333t  ♦  h22)  *  h16)) 
ob  =  osc(3)  +  oec(5)  +  osc(6)  -  dl  *  f(5)  *  (t8  -  It)  /  t6  -  .25t  *  fall  *  (6t  _ 

*  Hl7  *  (tl  ♦  2t  *  f(8))  -  h21  *  (tl  ♦  2t  ♦  2#  *  f(8)))  _ 

*  vlsl  *  h4  ♦  vla2  *  Hi  ♦  vl83  *  h5  ♦  agdg 

dl  =  dl  +  etaO  *  (vlel  *  h4  -  vll2  *  Hi  -  vle3  *  h5)  ♦  agdl 

e8q  =  de  *  de  ♦  dl  *  dl 

obc(4)  =  SQR(e8q) 

oac(3)  =  atan2(de  *  h19  ♦  dl  *  h18,  de  *  h18  -  dl  *  h19) 

osc(7)  =  It  -  2t  *  (di  *  di  ♦  dh  *  dh) 

fsi  =  SQR(lt  -  osc(7)  *  osc(7)) 

oac(6)  =  atan2(di  *  h8  ♦  dh  *  h7 ,  di  »  h7  -  dh  *  w8) 

oac(5)  =  08  -  osc(3)  -  osc(6) 
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ota2  =  1*  -  esq 
eta  =  SQR(«ta2) 

CALL  kepler(osc(3) ,  osc(4),  c3) 
osc(l)  =  C0S(c3) 
cl  =  If  -  obc(4)  *  osc(l) 
osc(l)  =  (osc(l)  -  osc(4))  /  cl 
osc(2)  =  eta  *  SIN(c3)  /  cl 

t84  =  If  ♦  08C(4)  *  osc(l) 

r  =  osc(8)  *  eta2  /  ts4 

w7  =  C0S(osc(6)) 

s8  =  SIN(o8c(6)) 

wl  =  C0S(obc(5)) 

s2  =  SIN(osc(5)) 

ubs  *  v2  *  obc(1)  «•  wl  *  osc(2) 

ubc  =  wl  *  osc(l)  -  w2  *  osc(2) 

u(l)  =  s7  *  ubc  -  b8  *  ubs  *  osc(7) 

u(2)  =  w8  *  ubc  ♦  w7  *  ubs  *  osc(7) 

u(3)  =  ub8  *  fsi 

ubc  *  osc(7) 
ubc  *  08c(7) 


v(l) 

= 

-w7 

* 

ubs  - 

v8 

v(2) 

= 

-w8 

» 

ubs  ♦ 

w7 

v(3) 

= 

ubc 

* 

fsi 

w(l) 

= 

u8  1 

'  1 

fsi 

w(2) 

= 

-w7 

* 

fsi 

w(3)  =  osc(7) 

ts2  =  eta  *  SQR(osc(8)) 

FOR  i%  =  1  TO  3 

vel(iy.)  =  (osc(4)  *  osc(2)  *  u(iy.)  +  ts4  *  v(iy.))  /  ts2 
NEXT 

vel(l)  =  vel(l)  +  cf3  *  r  *  u(2) 
vel(2)  =  vel(2)  -  cf3  *  r  *  u(l) 

'   Patch  for  PPT2  to  SATSTA  output 

elemu(l)  =  osc(8) ' semimajor  axis 

elemu(2)  =  osc(4) ' eccentricity 

'elemu(3)=  'time  of  last  perigee  passage 

elemu(4)  =  acos(osc(7)) ' inclination 

elemu(5)  =  osc(6)'node 

elemu(6)  =  osc(5)'arg  of  perigee 

elemu(8)  =  osc(3)'mean  anomaly 

maf  =  elemu(8) 

CALL  kepler(maf,  elem»(2) ,  eaf) 

sOf  =  SIN(eat) 

cOf  =  COS(eaf) 
'elemuf(8)  =  maf 
elemuf(3)  =  timef  -  maf  /  elemf(7) 

'  in-plane  position  and  velocity 

ecc2f  =  If  -  elemuf(2)  "  2 

sat.xf(l)  =  elemuf(l)  *  (cOf  -  elemf(2)) 

sat.xf(2)  =  elemuf(l)  *  SQR(ecc2f)  *  sOf 

sat.xf(3)  =  Of 

rf  =  elemuf(l)  *  (If  -  elemf(2)  *  cOf) 

edf  =  elemuf(l)  *  elemf(7)  /  rf 

sat.xdf(l)  =  -elemuf(l)  *  edf  *  sC» 

sat.xdf(2)  =  elemuf(l)  *  edf  *  SqR(ecc2f)  *  cOf 

sat.xdf(3)  =  Of 

'  equatorial  rectangular  coordinates 

CALL  euler(elemu(6)  ,  elemu(5)  ,  elem(4)  ,  sat.x(),  sat.xdQ) 

END  SELECT 
END  SUB 
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